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ABSTRACT 



The nutadonal stability of a dual-spin, quasi-rigid, axisymroetric spacecraft containing 
a driven rotor is analyzed. The purpose is to examine a revised energy-sink stability theory 
that properly accounts for the energy contribution of the motor. An inconsistency in the 
development disproves the existing energy-sink theory's assumption that the motor of the 
system contributes exactly enough energy to offset the frictional losses between the rotor 
and the platform. Using the concept of core energy, the revised stability criteria for a dual- 
spin, quasi-rigid, axisymmetric spacecraft containing a driven rotor is derived. An 
expression for nutation angle as a function of core energy over time is then determined. 
Numerical simulations are used to verify the revised energy-sink stability theory. The dual- 
spin, quasi-rigid, axisymmetric system presented by D. L. Mingori was chosen for the 
simulation. Equations for angular momentum and total energy were necessary to validate 
the numerical simulation and confirm aspects of the revised energy-sink stability theory. 
These equations are derived from the first principles of dynamics and are included in the 
analysis. An explicit relationship for core energy as a function of time does not exist. 
Various models postulating core energy are presented and analyzed. The numerical 
simulations of the computed nutation angles as a function of the postulated core energy 
compare well with the actual nutation angles of the system to confirm the revised energy- 
sink stability criteria. 
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I. INTRODUCTION 



A chronological review of early spacecraft types, and the stability criteria developed 
for them, provides a sufficient background for the fundamental concepts of this thesis. The 
revised energy-sink stability criterion is then presented, and an equation for nutation angle 
as a function of energy dissipation over time is developed. 

A. SINGLE SPIN SATELLITES 

1. Equations of Motion 

The earliest satellites took advantage of the fact that stability could be achieved via 
the 'gyroscopic stiffness’ of a spinning body. A preliminary dynamic model for the 
satellite could be achieved with Euler’s equations of motion for a rigid body. Eu 1 
moment equation can be written as 

M = N <U*- = B dlL + N a> B xh (lj 

dt dt 

In component form this becomes 

M\ — h\ + 0)2 /13 — 0)3 fi2 

M2 = fl2 + (D3 hi - 0)1 /13 (2) 

M3 = fl3 + CO\ /»2 — 0>2 h\ 

Simplification of the model is accomplished by assuming that it is axisymmetric, that the 
body-fixed axes coincide with the principal axes (correct for simple spacecraft with 
/ 1 =/2)> and that the body is in a torque-free environment (a valid approximation used 
throughout this thesis). The spacecraft is then represented by Figure ( 1 ) and the equations 



1 



of motion are 



0 = 1 1 (0\ + (/3 - / i ) a>2 CO3 

0 = / 2®2 + (/l- 73) CO 3 CO 1 ( 3 ) 

0 = I 2 CO 3 



b 3 




Axisymmetric Rigid Body 
Figure ( 1 ) 



The angular velocity vector, angular momentum vector, and the kinetic energy may be 
expressed respectively as 



N-B 

o) = a>i bi + 0)2 b2 + C03 b3 


( 4 ) 


h = I\ coi bi + 12 C02 b2 + h C03 b3 


( 5 ) 


E = 1 (/1 ( 0 \ + h CO 2 2 + h CO 3 2 ) 


( 6 ) 



A simplification of the notation can be made. Let I\ = h = A and h-I s • From the third 
line of Equation ( 3 ) it can be seen that the angular velocity about the spin axis b3 is 
constant, therefore 03 = co s . The transverse angular velocity components interchange 
between the bi and the b2 axes but the magnitude is constant, so one can let 



2 



co, = co i bi + 002 t>2- Therefore, 



'~B 

0) = CO, + CQ S b3 


a) 


h = /, c5, + / 5 q) s b3 


(8) 


h 2 = I t 2 co 2 + I 2 co 2 


(9) 


E = !•(/, a>, 2 + /,<»/) 


(10) 



Because the motion is torque-free, |h| is constant and h is fixed in space. Because there are 
no energy sources or sinks, E is also constant. Finally, the above conditions result in 

co^l being a constant. 

An expression for the nutation angle may now be developed. The orientation of 
the body axes with re: ~t to an inertial reference frame is desired to provide a measure of 

the body's dynamic behavior. Because the angular momentum vector is fixed in space, the 
nutation angle is defined as the angle between the body-fixed axis about which spin is 
desired and the angular momentum vector, and can be expressed in one of the following 
forms 




In these equations the first and second terms correspond to the general case with spin about 
the b 3 axis, and the third term uses the simplified notation to describe an axisymmetric 
body. In the special case of spin about only one principal axis of an axisymmetric body, 
the angular momentum vector and the angular velocity vector will lie on the spin axis. With 
spin components on two or more principal axes, the three vectors will not be coincident, 
although they still will lie on the same plane. 
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2. Single Spin Satellite Stability 

A torque-free, axisymmetric rigid body with the body axes coinciding with the 

principle moments of inertia will be stable about the axis of either the maximum moment of 

inertia or the minimum moment of inertia. To prove this, one begins with an arbitrary rigid 

body. The body is given the initial condition of steady angular velocity, ct)o, about a 

principal axis, and is then perturbed slightly. It is assumed that the angular velocities about 

the other axes are small, and are approximately the same order of magnitude as the 

perturbation (o)i = 0) 2 = e). The system will be considered stable if the perturbation does 

not increase over time. Given an initial angular velocity with a perturbation, 
N-+B 

0) = o>i bi + 002 b2 + (0)0 + e) b3, and given arbitrary inertias / 1, I2, I3, Euler's 
equations of motion can be written as 

0 = /1 ©1 + (I3 - h) o ) 2 {coo + e) 

0 = / 2 0>2 + (/1 — 13) Ct)\ (Q)q + £) ( 14 ) 

0 = 13 a >3 +(/ 2 -/ i ) g>i (02 

The equations are linearized by neglecting terms of magnitude e 2 . Rewriting the equations 
by eliminating the terms e coi, e o> 2 , and o)\ n > 2 results in 



72-/3 

G>1 = , - (02 0) 0 

h 


( 15 ) 


(/3-/1) 

0>2 = ■— r — 0>i 0) 0 

h 


( 16 ) 


(O3 = 0)q + £ = £ = 0 


( 17 ) 



From Equation ( 17 ), one can conclude that e= constant. By differentiating Equation ( 15 ), 
and using Equation ( 16 ) to eliminate 6) 2, one gets 

< “ l+ ( (,3 ~ / / 2 | ) l' 3 — M o 2 )^. = 0 (18) 

Similarly, from differentiating Equation ( 16 ), and using Equation ( 15 ) to eliminate 6) 1, one 
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arrives at 



(i>2 + — ^O)o 2 ja)2 = 0 (19) 

These equations are identified as second order, linear, ordinary differential equations with 
constant coefficients. The general solution for these differential equations are 

e‘>' + e-'>' (20) 

o >2 = K 2 c i r t + K 4 c~ i r t (21) 

where 

r = Z2^r ( 22 ) 

If > is imaginary, coi and ©2 will increase without bound over time, and the motion will be 
unstable. Stability is achieved if y is real. The first case occurs when the maximum 
moment of inertia is about the spin (b3) axis; then I 2 > I\ , h > h, and (I 2 - I2) {h - A) > 0. 
In this case the inertia ratio I s / /*, defined as the inertia about the spin axis over the inertia 
about a transverse axis, is greater than one. The second case occurs when the minimum 
moment of inertia is about the spin (63) axis; then I3 < I\ , 1 3 < 1 2, and (1 3 - h)(h ~I\)>0 
as well. Here the inertia ratio is less than one. In both of the above cases y is real and the 
motion is stable. However, if I3 is the intermediate moment of inertia about the spin (b3) 
axis, then (73 - 12) {I3 - 7 i) < 0, y is imaginary, and the motion is unstable. 

The previous model cannot be applied to a satellite since the assumption of a rigid 
body cannot be extended to the spacecraft. Structural elasticity, liquid propellant slosh, 
etc., cause energy dissipation in an actual spacecraft. This spacecraft can be generalized by 
a quasi-rigid body with an unspecified energy damper mechanism. A priori, one can 
conclude that energy in the above system will dissipate until the minimum energy state is 
reached. The kinetic energy of the system with spin about the principal axis with maximum 
moment of inertia and spin about the principal axis with minimum moment of inertia can be 
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written respectively as 



E = 1-AL 

21 max 



£ = 



- 1 hL 



21 



I m/TT about b 3 
/min about b 3 



(23) 



mm 



The angular momentum is constant in the torque-free case. Thus the minimum kinetic 
energy state occurs when the rotation is about the axis of the maximum moment of inertia. 
Therefore, a quasi-rigid body is stable only when it is spinning about its major axis, with a 
corresponding inertia ratio that is greater than one. 

A relationship can be established between the time rate of change of the nutation 
angle and the energy dissipation of a quasi-rigid, axisymmetric body. One must assume 
that the angular momentum and moments of inertia of the quasi-rigid body do not change 
appreciably from a comparable rigid body. For the generalized model, with arbitrary 
inertias I\=h, and / 3 , Equations (5) and (6) are substituted into Equation (12) to obtain 



sin 2 6 = 



\ (2/ 3 E-/t 2 ) 

h 2 



(24) 



\h-h) 

The time rate of change of the nutation angle is determined by taking the derivative of the 
above equation 

1 2/i/a 



6 = 



' Ed total 



(25) 



sin ( 20 ) (/ 3 ~l\)h 2 

where the only rate of change of energy E is attributed to the damping mechanism and is 
written as Ed total to emphasize this point. Because Ed total is negative, the nutation angle 
will decrease only if / 3 is greater than I\. This reaffirms the previous conclusion that a 
quasi-rigid body is spin stabilized only about the axis of maximum moment of inertia. The 
foregoing development is referred to as the energy-sink method. 
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B. DUAL SPIN SATELLITES 

The logical progression from the single spin satellite was to incorporate a de-spun 
platform. This permitted the replacement of the low gain omnidirectional antenna with 
directional antennas for communication satellites, and a more capable spacecraft for 
scientific observation. A simple control system about the b 3 axis would maintain the 
platform rotating at a constant relative rate with respect to the rotor ( and would usually 
have the platform rotate at the earth’s rotational rate). Initially, the platforms were 
sufficiently small, and the overall dimensions of the satellite were such that the inertia ratio 
would be greater than one. For this type of satellite, the previously developed theory 
proved adequate. However, as satellites continued to grow in size, the launch vehicle 
shroud diameter became a constraint. In order to provide the size spacecraft needed to 
satisfy mission requirements and still fit within the shroud, a spacecraft with an inertia ratio 
of less than one {I s total! It total < 1) would need to be built. From the previously 
developed theory, a spacecraft with an inertia ratio of less than one was believed to be 
inherently unstable. It was not until the development of the energy-sink theory for dual- 
spin, quasi-rigid, axisymmetric spacecraft containing a driven rotor, developed 
simultaneously by V. D. Landon (unpublished work) and A. J. Iorillo [Ref. 1], that a 
spacecraft with an inertia ratio less than one was considered feasible. Several rigorous 
stability analyses using the equations of motion for specific dual spinners have been 
performed by P. W. Likins [Ref. 2], D. L. Mingori [Ref. 3], and others, to validate the 
energy-sink theory. The difficulty of a rigorous analysis is in accurately modeling all the 
forms of energy dissipation. A more general and practical approach was required to 
determine stability, and the energy-sink theory proved suitable. The development of this 
theory is as follows. 
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The spacecraft, shown in Figure (2), is assumed to fulfill the following conditions: 

• Both the rotor and the platform are axisymmetric 

• Both the rotor and the platform are quasi-rigid 

• The damping mechanisms do not significantly alter the energy value, although 
the mechanisms will contribute to the energy rate 

• No external torques are applied 

• The only relative motion is spin about the b 3 axis 

• The motor, driven by the control system, inputs just enough energy to exactly 
offset the shaft frictional losses, maintaining a constant relative angular 
velocity between the rotor and the platform 




Dual-Spin Quasi-Rigid Axisymmetric Spacecraft 
Figure (2) 
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The magnitude of the angular momentum and the kinetic energy of the dual-spin system can 
be expressed respectively by the following equations 

h 2 = I? total CO , 2 + (is (O s + Is (Osf = I? total <*>} + 0 + h total (O/f (26) 

E = ~ It total (O 2 1 s 0) 2 I s (O s = ^It total (0^ Is total (Os + 2 ’/jCT^+/ j <U j 0 (27) 

where 0 is the relative rotation rate of the rotor with respect to the platform. Although the 
equations actually represent a rigid body system, they are also applicable to the quasi-rigid 
system because of the above assumptions. If the damping mechanisms do make significant 
contributions to the energy of the system, then Equations (26) and (27) do not hold, and 
the energy-sink criterion will not apply. Additionally, the potential energy of the system 
(for example the energy stored in the spring of a mass-spring-dashpot damper) will not 
make a significant contribution to the total energy of the system. Therefore, for the 
remainder of the thesis, it is assumed that the kinetic energy of the system is effectively the 
total energy of the system (e = E tota i). Because the system experiences no external torques, 
angular momentum is conserved. Because the motor contributes no energy to the system, 
the time rate of change of the kinetic energy E is represented by only the quasi-rigid energy 
dissipation Ed total » and is negative. Differentiating the above two equations with respect to 
time, one obtains 

o = rftotal tot (Ot + (is (Os + I s' (Os') (is (0 S + I s' (Os) (28) 

E — Eq total ~ It total (Ot (Ot + Is (Os (Os Is (Os (Os (29) 

Eliminating the common term (O t cb, by combining the above equations 

(Os + IJ_ (Os') (is (Qs + IJ_ (Qs 
It total 

One may now define the inertial nutation frequency, Xq, the rotor nutation frequency. 



+ I s (0 s (b s + 1 s ' (Os' (bs' (30) 
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A, and the platform nutation frequency, A\ as follows 



* 7j 0) s + I s G) s 

AO = - 

h total 



( 31 ) 



-3 o (Is -It total) CDs + Is &s 

A = Aq- CO s = 

*t total 



' — ~ ll total) 

It total 



X = Ao - fl>/ = 

The nutation angle for a dual-spin system is defined as 






By imposing the condition 



(32) 

(33) 



(34) 



Xo > 0 (35) 

the analysis of nutational motion is restricted to the following region without any loss of 
generality 

0<*<f (36) 

Incorporating the nutation frequency terms, the equation for energy dissipation is written as 

E = Ed total = Ed + Ed = -I s Xco s - 1 / X' co s ' (37) 

Recalling the assumption that the rotor and platform are uncoupled about the b 3 axis, we 
may incorporate Equation (37) into the reaction torques which tend to change the angular 
rates 




// (Os' 




X' 



(38) 
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Combining Equation (30), (37), and (38), one arrives at the transverse rate equation 



h total <*>t G>t = Ao (— + 

. r 



( 39 ) 



Differentiating Equation (12) and substituting Equation (39) into it, the time rate of change 
of nutation angle as a function of energy dissipation rates is 

8 = 2 I, touU Ao (go (40) 
sin (2 d)h 2 \X A '/ 

The energy-sink equation for a single spin-stabilized body is obtained by letting A and A' 
become Ao and Ed + Ed become Ed total for a single body. The definition of stability 

requires the nutation angle 8 to remain constant or decrease as a function of time, so that 8 
is zero or negative. Because the factors outside the parenthesis on the right hand side of 
Equation (40) are positive, the stability criterion for a dual-spin, quasi-rigid, axisymmetric 
system becomes 




One of the following cases will guarantee stability 



(41) 



1) A>0 and A'>0 



A > 0 , A' < 0 , and 


Ed’ 


< 


Ed 




A' 




A 


A < 0 , A' > 0 , and 


Ed 


> 


Ed 




A' 




A 



A specific example would be the model of a typical communications satellite, a prolate 
dual-spinner possessing an inertia ratio of less than one. In general, the rotor nutation 
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frequency, expressed as 

X — — co — \LlZ-Ll {£f£fl jj. (42) 

It total It total 

would be negative because (I s - It total) is negative and co s » co/ if co/ is rotating at the 
earth's rotation rate. Thus, energy damping in the rotor, Ed, would be destabilizing while 
energy damping in the platform, Ed', would be stabilizing. It is from this result that 
satellites will have a damping mechanism placed on the platform to improve nutational 
stability. Such a damper is called a nutation damper. 
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n. PROBLEM DEFINITION 



A. OVERVIEW 

The existing energy-sink theory relies on several assumptions, perhaps the most 
important relating to the driven rotor. As previously stated, it has been assumed that the 
motor, driven by the control system, inputs just enough energy to exactly offset the shaft 
frictional losses, maintaining a constant relative angular velocity between the rotor and the 
platform. In actual systems, contrary to this assumption, the motor may add or remove 
energy from the system, depending on the dynamics of the spacecraft. Consequently, a 
revised energy sink stability theory, properly accounting for the energy contribution of the 
motor, is derived. The revised theory, based on the concept of core energy, will remain 
consistent with the existing energy-sink stability criterion. Continuing, an equation for 
nutation angle over time, as a function of core energy, is developed. Given a postulated 
energy dissipation function modeling the nutation dampers, structural elasticity, fuel slosh, 
etc., one can accurately predict the nutation angle behavior. Numerical simulation of D. L. 
Mingori's dual-spin, quasi-rigid, axisymmetric system containing a driven rotor [Ref. 3] is 
used to validate the revised energy-sink stability theory. The predicted nutation angle, 
based on this revised energy-sink theory, and the postulated energy dissipation function, is 
compared to the exact nutation angle of the Mingori system. By using a suitable postulated 
energy dissipation function, one can achieve excellent agreement between the predicted and 
the exact nutation angle. I. Michael Ross [Ref. 4] performed this analysis on a dual spin 
system with a damper on the platform only. The remainder of this thesis will use the same 
analysis, but on the Mingori system with dampers on both the rotor and the platform. 
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B. INCONSISTENCY OF THE ENERGY SINK THEORY 



There exists a contradiction between the existing energy-sink theory and the nutation 
angle derived from it. This will provide the motivation for developing a revised energy- 
sink theory and an alternative equation for nutation angle over time. From before, the 
existing energy-sink stability criterion can be expressed as 



^> + ^<0 
A A' 



(41) 



and the nutation angle for the dual-spin system was defined as 



0 Jh-b 3 \ ,lh3 + h 3 '\ , (l s co s + I/ co s 

6 - cos -1 — in — = cos -1 — n — = cos -1 m 



N 



|h| 



M 



(34) 



As previously stated, for a prolate dual-spinner ( I s , ota i / /, tot ai < 1 ), energy dissipation in 
the platform is stabilizing and energy dissipation in the rotor is destabilizing. The angular 
velocity of the platform about the spin axis, (o/, can be expressed in terms of the angular 
momentum and the kinetic energy of the system. Combining Equation (26) and Equation 
(27), one arrives at 



0 )/ = - 



h_Q 

Is total 



+ - / 1 Is CT — 2. It total E + I s {It "total ~ 7s) 

> Uj total ' I s total Vt total ~ Is total) 



Substituting this expression into Equation (34) results in 

6 = cos -1 (± Vq) 



(42) 



(43) 



where Q is represented as 



( 2 E-I s o 2 ) I s total + Ij a 1 


Is total 


^^Js total) 


It total ” Is total 


' h total i 





(44) 



Initial conditions at t = 0 will determine the correct sign, with continuity considerations 
maintaining the sign for all of t > 0. Additionally, no external torques are applied to this 
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system, resulting in |h| being constant Differentiation of Equation (43) results in 



6 = - 



1 d (+ 4q\ = ± — 1 1 Ed total) Is 



total 



Vi-(±Vs) 2 ‘* 



sin 6 4 q ^[i -(tiflatj 

V It total ' 



( 45 ) 



From the definition of nutation angle. Equation (34), and the condition imposed on it. 
Equation (36), the positive sign must be chosen in Equations (43) and (45). An important 
observation is made at this time. Choosing the positive sign will result in a positive rate of 
change in the nutation angle, indicating an unstable condition. The relative rotor spin rate is 
an independent variable, and is arbitrarily selected here as a constant value over time. 
Therefore, energy dissipation in a prolate dual-spin spacecraft will make the nutation angle 
increase, regardless of whether the dissipation is in the platform or in the rotor. This is not 
consistent with the stability criterion of Equation (41). Thus, the existing energy-sink 
stability criterion contradicts itself. 



C. CORE ENERGY AND ENERGY DISSIPATION 

The existing energy-sink stability criterion does not properly account for any energy 
that may be provided by, or absorbed by, the motor. To accurately represent the system, 
the total rate of change of energy must be written as 

E = E D total + W (46) 

where W is the rate of work due to the motor, and may be either positive or negative and 
the rate of change of energy due to dissipative elements can occur in either the platform or 
the rotor. Recall that the kinetic energy of the dual-spin system was expressed in Equation 
(27). If the work due to the motor torque as a function of time is written in analytical form, 
then the time rate of change of the energy of the system due to all dissipative elements, 

Ed total » can then be expressed solely in terms of the quasi-rigid parameters of the dual-spin 
system. With this expression, the condition that Ed total ^ 0 will result in the required 
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stability criterion. The difficulty arises in that to get the work due to the motor torque W 
(or W), one needs to know the exact dynamics of the dissipative mechanism. In the 
development of the modified energy-sink stability. Equation (46) will be used to determine 
the expression for Eq tola [, and ultimately derive the revised stability criterion and nutation 
angle equation. 

Additionally, the existing energy-sink theory can be shown to be incorrect for both 
the case of total energy decreasing and for the case of total energy increasing. For 
example, if the rotor is rigid and total energy decreases, Equation (41) predicts that the 
system will be stable, but Equation (45) predicts that it will be unstable. Allowing the total 
energy to increase would reverse the conditions, but still show a contradiction between 
Equation (41) and Equation (45). 

A modification of the energy-sink stability theory and the associated expression for 
nutation angle is now presented. The development of the theory is from I. Michael Ross' 
unpublished notes. The basis of the new theory is centered on the core energy of the 
system. As defined by Hubert [Ref. 5] 

Core energy is the total energy of the spacecraft minus the portion of the rotor 
energy that is due to the relative rotation between the rotor and the platform. It 
is assumed that the mass, inertia, and motion of the damping device are 
sufficiently small that its energy is negligible relative to the spacecraft core 
energy. The damper will be treated as an undefined 'energy sink' for the 
purposes of the energy sink analysis. 



From the above statement one can define the core body as that body whose inertial 
dynamics are selected for analysis. Hubert defines the platform as the core body. The core 
energy is simply the rotational kinetic energy of a fictitious rigid body that possesses the 
inertia properties of the entire dual-spinner but moves in inertial space exactly like the core 



16 



body. For a dual-spin spacecraft the platform core energy is defined as 

Ec = total & i + 2 ^ 2 total 0)2 + 2 toUd 0)3,2 = 2 ^ totai 0)2 + 2 I* tolal 0)5 2 
The center expression is for an arbitrary dual-spin spacecraft with the spin axis about the b 3 
axis, and the right expression is for a dual-spin, axisymmetric spacecraft with the 
simplified notation. Extending this concept to the rotor, the rotor core energy is defined as 

E C — 2 " ^1 tota l ® 1 2 " ^2 total ^2 ^ ^3 total ^ 3 ^ ~ toia2 + 2 tota ^ (^ 8 ) 

Necessary to the development of the modified theory is what will be termed the 
Separation Axiom. This is when analysis is first performed with the rotor considered rigid 
and the platform quasi-rigid. Euler's equations are written for the rotor, and through 
manipulation, an equation is derived relating the torque on the quasi-rigid platform solely in 
terms of platform variables. Then the platform is considered rigid and the rotor quasi-rigid. 
An equation is derived relating the torque on the quasi-rigid rotor solely in terms of rotor 
variables. These two separate equations are then combined and applied to a system in 
which the rotor and the platform may both be quasi-rigid. 

The case of a rigid rotor with a quasi-rigid platform is first analyzed. Because the 
rotor is rigid, the torque applied by the motor to the rotor is the net torque on the rotor and 
is determined from Euler's moment equations. For the case of the axisymmetric rotor, 

Tr - Trim - I s o) s (49) 

One can observe that Tp/m = - Tr/m , but 7> * Tr since the damping mechanism 
contributes additional torques to the quasi-rigid platform. The rate of work needed by the 
motor torque to maintain constant relative motion between the rotor and the platform is 

W = Tr g = I s G(O s = I s <7(6)/ + o) (50) 

By substituting the platform core energy, Equation (47), into the kinetic energy expression. 
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Equation (27), kinetic energy may be expressed as 

E = Ec' + Li s g 2 + I s O)s' o (51) 

The above equation is differentiated to arrive at the time rate of change of the kinetic energy 
of the dual-spin system 

E = Ec + l s o o+ I s (era)/ + (Os o) = Ec + I s o{co s ' + oj + I s 6 (52) 

Substituting into this equation the rate of work done by the motor torque, Equation (50), 
one gets 

E = E c ' + W + I s (Os' 6 (53) 

Comparing this to Equation (46), it can be seen that 



Eq total — Ec + Is co s o (54) 

Taking the derivative of Equation (47) to get the time rate of change of the platform core 
energy 



Ec = It total CO, 6), + Is total CO s ' CO s ' (55) 

and then substituting this into Equation (54), one arrives at the total energy dissipation of a 
rigid rotor, quasi-rigid platform system 



Ed total — It total CO, CO, + I s total CO s (O s + I s CO s <7 

= It total CO, 0), + I s CO s ' COs + V CO, COs + h CO s 6 



(56) 



Because the system has no external forces applied, it remains torque-ffee. Thus, Equation 
(28) can be used to eliminate the co, co, term and arrive at 



Ed total — {is C0 S + Is COs ) 



{I t total — Is total) CO s Is CS 
It total 



(57) 
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Noting that the platform nutational frequency. Equation (33), can be written as 



A' = Ao -©/ = 



Is (O s + I s 0 ) s \ , _ {Is total - total) &s + Is & 

I “ W s — 



It total 



It total 

then Equation (57) can be written as 

Ed total = (is (Os + Is (Os') (- A') = Ao I,(- A') 
Referring to Equation (34), the nutation angle can be written as 



(58) 



(59) 



cos 6 = 



(Is (Os +1/(0/ 

l N J 



(60) 



Taking the derivative and comparing it to the rate of total energy dissipation. Equation (59), 
the following relationship can be established 



- 0 sin B = 



I s 6) s + I/6)/\ -E D tota i 



\ l h l I A'|h| 

Because the rotor is rigid, all energy dissipation will occur in the platform, such that 

En 

— = |h| 6 sin 6 
A' 



(61) 



(62) 



Equations (61) and (62) can be rewritten by including the motor torque. Equation (49), as 



|h| 6 sin 6 = T r + I/6)/ = 

X 



(63) 



Because Tpim = - Trim (action - reaction pair), the final relationship is written as 

Trim = — + 1/ (o/ (64) 

A' 

This equation describes the motor torque on the quasi-rigid platform as a function of 
platform variables only. It can be seen at this point that the classical analysis can be 
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achieved by assuming no torque is applied by the motor, resulting in Tr/m = 0 and 

— = -h' co/ (65) 

X' 

This analysis can now be performed on the system containing a quasi-rigid rotor and 
a rigid platform. It is observed that the selection of the body to be the rotor and the body to 
be the platform is completely arbitrary, and no physical distinction exists between the two 
bodies. Therefore, by analogy, the equation describing the torque on the quasi-rigid rotor 
as a function of rotor variables must be 

Tr/m = — + I S 0) S (66) 

X' 

Now let Tp/M and Trim represent the motor torques on the platform and on the rotor 
respectively for the system containing both a quasi-rigid rotor and a quasi-rigid platform. 
Then the separation axiom would require the following two conditions 

Tp,M = Tp/M = — + // (Os (67) 

X' 



Tr/m ~ Tr/m - — + h 

X 

Once again, since the system is an action - reaction pair 



Tp/M + Tr/m - 0 



and then 



Ed | Tp 
X X' 



= ~(l s Q) s + 1 s' CDs') = |h| e sin e 



(68) 



(69) 



(70) 
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The stability criterion dictates that 6 remain zero or be negative, thus 



^+— £ 0 0 
X X 



This is seen to be the existing energy sink criterion. Therefore, despite the presence of 
motor torque, the existing stability criterion still applies. Referring to Equations (67) and 
(68), the second term in the equations would represent the torque applied by the damping 
mechanism to the platform and the rotor respectively. Equations (67) and (68) state that the 
motor torque minus the damper torque will equal the net torque of the platform and the 
rotor respectively. 

In determining the revised energy-sink stability theory, the energy dissipation 
equation for the system with no energy contribution from the motor. Equation (37), must 
be rewritten to account for the motor torque. Thus 



E — Eq total W — — /$ A Ct)j - /j 1 COj 
Substituting in Equations (67) and (68), one arrives at 



(71) 



£ = t: Dlolal +w = Ap-rjJ+r 



= e d *e d ’-(xt' riu *xt' pim ) 

By assigning the motor torque values as 



U' 



P/M I 



(72) 



t m - T r/m T pim (73) 

then the above equation can be rewritten as 

E = E D total + W = E D + E D ' + T M (X'-l) (74) 

Because the total time rate of change of energy dissipation equals the rate of change of 
energy dissipation in the rotor plus that in the platform, the rate of work due to the motor 
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torque must be 



W = T M (X'-X) = T m o 



(75) 



Referring to the example from Chapter 1, the prolate dual-spin communications 
satellite, typically the rotor nutation frequency would be negative and the platform nutation 
frequency would be positive. For rotor spin-up Tm> 0 and the motor torque would add 
energy to the system, and for rotor spin-down Tm< 0 and the motor torque would remove 
energy from the system. For the arbitrary system, the motor torque, Tm , and the sign of 
the term X' -X (= o), would determine whether the motor adds or removes energy from 
the system. 

The rates of change of the energies of the system may now be represented. Rewriting 
Equation (46) to determine the total energy dissipation of the system 



The rate of work due to the motor torque can be expressed by combining Equations (75) 
and (68) (or Equation (67) ) to arrive at 



The rate of change of the kinetic energy of the system as a function of platform core energy 
and system parameters was determined previously as 



Substituting the above two equations into Equation (76), the expression for total energy 
dissipation of the system may now be written as 



Ed total — E — W 



(76) 




(77) 



E = E c ' + i s o (©/ + <y) + i s 



(52) 




( 78 ) 
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Because 6) s ' + a= cb s , this equation can be further simplified 



Ed total = E C ' + Is (0/ 6- — o (79) 

A 

Comparing this equation, representing the system with a quasi-rigid rotor and a quasi-rigid 
platform, with Equation (54), representing the system with a rigid rotor and a quasi-rigid 

platform, reveals an additional term, - ^ o. The first term of Equation (79) represents the 

A 

rate of change of core energy, with the platform as the core body. The second term will 
account for the change in energy associated with a change in the relative rotation rate of the 
rotor with respect to the platform. The final term accounts for the energy loss due to the 
quasi-rigidity of the rotor that is not represented in the platform core body expression. 

The case that will be analyzed is that which occurs when there is a constant relative 
rotation rate between the platform and the rotor. For the remainder of the analysis, let 

<7=0 (80) 

and the total rate of energy dissipation becomes 



Ed, cat = Ec ’-=22- (81) 

X 

Further simplification can be achieved by noting that a = A' - A and from Equation (37) 
that Ed total = Ed + Ep'. Therefore 

E d + E d ' = E c ' - E D (82) 



which reduces to 



Ed + Ed_ = Ec_ 
A A' A' 



(83) 



A similar development can be performed using rotor core energy vice platform core energy. 
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The result is 



Ed + Ed_ = Ec 
X X' X 



( 84 ) 



There is an expected symmetry between Equations (83) and (84), due to the arbitrary 
assignment of one body of the system as the rotor, and the other body as the platform. To 
confirm the results of Equations (83) and (84), one must prove that when o is constant 

^ — (85) 

X X' 



Taking the derivative of the rotor core energy, Equation (48) 

Ec = It total ®t Mt + Is total ®s 

= It total O), 6), + I s total ( CO / + o) (ft)/ + c) 



and similarly, for the platform core energy 

Ec = It total ®t Is total (87) 

Substituting Equation (87) into Equation (86) and noting that <7=0, 

Ec — Ec + Is total G (88) 

Multiplying through by the platform nutation frequency 

E c X' = E c ' X' + I s total % o ft)/ (89) 

and recalling that X' = X + o, 

Ec X' — Ec X + Ec <J + Is total % Q &>/ (90) 

which is the same expression as Equation (85) if it can be proven that 

Ec + Is total % ft)/ = 0 (91) 



Eliminating the transverse angular velocity of Equation (47) by substituting in Equation 
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(26), one arrives at 



£c ' = 2 



U 2 - 



(/j Q) s + Is (Os Y | 1 , 

I, total I 2 5 



total (Os 



,2 



(92) 



Taking the derivative 

Ec 



fr + V «/) + V «/) . , „ / • / // ***/ 

-C — } *" *s total W s CO s 

*t total *t total 



(93) 



and simplifying 

2 

5 2 COs Ois + I_s I s' (Os &>/ + /j // (Os V <*>/<»/ 

It total 

Noting that / 5 total - I s + Is and recalling that a)* = co s ' because <x = 0, the above 
equation can be reduced to 




) + Is total I t total (Os (Os 



(94) 



E C — ~ Is total 



(/j (Os + Is (Os It total (Os ) (Os 



It total 



(95) 



finally, invoking the definition of the platform nutation frequency, one arrives at 



E C — ~ Is total A (O s 



(96) 



Therefore, Equation (91) holds and the revised stability criterion can be expressed as 



Ec _ Ec_ _ Ed ( Ed < Q 
A A' A A' 



(97) 



A few remarks can be made concerning this stability criterion. The third expression is the 
existing energy sink stability criterion. This criterion must equal the stability criterion as a 
function of the rotor core energy which must equal the stability criterion as a function of the 
platform core energy. It can be seen that one no longer needs to know the energy 
dissipation rates in the platform and in the rotor to determine stability. By knowing or 
postulating the core energy over time, the stability of the system can be determined. 
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Continuing with the prolate dual-spinner example of Chapter 1, any one of the above 
expressions apply. For a stable system, the rotor core energy will be positive, and increase 
over time. Additionally, the rotor nutation frequency will be negative, resulting in a 
negative expression for the rotor core energy stability criterion. The platform core energy 
will be positive, but will decrease over time. The platform nutation frequency will be 
positive, resulting in a negative expression for the platform stability criterion. According to 
Equation (97), the platform and rotor stability criteria must equal one another. From the 
numerical simulation of a dual-spin quasi-rigid axisymmetric system, the rotor and platform 
core energies as a function of time will be determined and graphed. 

D. NUTATIONAL MOTION 

The development of a modified expression for the nutation angle as a function of time 
may now be presented. The actual nutation angle of the system is defined as 6. The 
nutation angle as a function of platform core energy will be represented by rf, and the 
nutation angle as a function of rotor core energy will be represented by rj. The derivation is 
similar to the one previously given in this chapter, except that the total energy of the system 
has been replaced by the core energy of the system to eliminate the transverse angular 
velocity, to, . The derivation for the platform core energy will be shown. From before, the 
angular momentum of the dual-spin, quasi-rigid, axisymmetric system is 

h 1 = I? total °>t 2 + (h CDs + Is CO s ') 2 = I? total CO t 2 + (l s O+ I s total 0)/Y (26) 

Combining this with the platform core energy, and solving for the platform angular velocity 
about the spin axis 

I s total (I s total ~ 1 1 total) G ) 3 — 2 7 s total I s G +I 2 0’ 2 + 2Ec It total ~ h 2 = 0 (98) 
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and 



«*' = Tj ^ r±Vg( W 1—1— 

Vl to/a/ “ -*5 to/a// \*/ to/a/ “ *5 total) *s total 



(99) 



where 



£?' = 2 £ c ' /j total ( 1 - _ /,2 j f ] _ ZiJ2£ai| + \Lj 2 isl\ j 2 <7 2 (100) 

V 1 1 total ) \'t total ) V *t total ' '*t total ' 

The initial conditions at t - 0 will determine the proper sign of Equation (99). As before, 
continuity will ensure this sign for all t > 0. Substituting Equation (99) into Equation (11) 



A 1 * total (Oi+l s o\ , 

Tj = COS 1 I — J = COS 1 



It total 



L[ho±-Jfr\\ ( 101 ) 



J t total ~ I s total. 

The time rate of change of the nutation angle will determine the stability. Differentiating the 
above equation results in 



W — !?£. J. ( ll total \ 

' sinri h [iVe 7 ] 



(102) 



A stable system would require that the lower sign be chosen for the radical, thus 



77 = cos 



-1 



It total 



\\i s g-*Jq 



\ ft total ~ I s total, h 

When this sign selection is applied to Equation (52), one gets 




(103) 



(A total - Is total) + /, a = - VcF [LuotaA 

Oj total I 



(104) 



This can be rewritten as the well established dual-spin stability condition, as written by P. 
C. Hughes [Ref 6], 

(t t total ~ 1 s total) ^3 + 1$ & — 0 (105) 

It is important to note that this analysis does not produce a contradiction to Equation (41). 



27 



In a similar manner, the modified nutation angle can be derived with respect to the rotor as 



n = cos-»([ , '-“*f — •ljL[v<r-V5]| 

\l/f total — *s total J ft I 



(106) 



where 



Q = 2 E c I stota i[\ -^-i2£fll|_^ 2 (£j_£2£fllj(l + q2 (107 ) 

' ’t total' \*t total' ' 't total' v/f total' 

Therefore, if one is given the core energy or the postulated core energy as a function of 
time, the nutational motion can then be determined. This leads to an extremely important 
conclusion. By determining a sufficiently accurate postulation of the core energy of a dual- 
spin, quasi-rigid, axisymmetric spacecraft over time, one can predict the nutational 
behavior and the stability of that spacecraft. Numerical simulation of such a system will 
confirm this conclusion. 
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HI. DYNAMICAL EQUATIONS 



A. MINGORI DUAL-SPIN SYSTEM 

A specific model is required to validate the modified energy-sink stability theory. The 
development of the previous chapters was completely general in nature and applies to any 
dual-spin spacecraft with an arbitrary damping mechanism. D. L. Mingori’s dual-spin 
system [Ref. 3 ] provided the needed model required to validate the proposed stability 
theory. Figure ( 3 ). Additionally, the complete non-linear equations of motion were 
presented by Mingori; however, the expressions for the important parameters in attitude 
dynamics, angular momentum and kinetic and total energy, were not presented in his paper 
and had to be derived before any analysis could be performed. 

The Mingori system is comprised of two symmetric rigid bodies, the lower which 
shall be defined as the rotor, and the upper body shall be defined as the platform. By 
convention, all terms referring to the platform will be the same notation as that of the rotor, 
except that they will carry the prime mark. Both the rotor and the platform have coordinate 
axes fixed to the body and located at the respective centers of mass. The distance between 
the centers of mass is specified by L. The entire spacecraft has coordinate axes fixed to 
the spacecraft center of mass, denoted by the double prime coordinate axes, and rotating in 
the same manner as the rotor coordinate axes. The coordinate axes b3, b3', and b3"are all 
collinear. The spacecraft center of mass will vary along the b3 axis as the point masses in 
the rotor and the platform oscillate. The only relative motion of the platform with respect to 
the rotor is angular rotation about the b3 axis. A motor driven by a control system 
maintains a constant relative rotation rate o. The angle between a line parallel to bi and 
passing through the platform center of mass, and *>r is represented by y = ot . 
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Mingori Dual-Spin Quasi-Rigid Axisymmetric System 
Figure (3) 
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Each body contains a mass-spring-dashpot mechanism. An actual spacecraft 
undergoes damping from various mechanisms. Undesired damping can occur due to 
structural elasticity and liquid propellant slosh. Nutation dampers are incorporated into 
spacecraft to improve the nutational motion. The dual-spin axisymmetric system with 
mass-spring-dashpot dampers cannot accurately model an actual spacecraft, but is used to 
illustrate and validate the theory. A description of the rotor is as follows. The mass- 
spring-dashpot mechanism can be modelled by a particle mass m\ attached to a spring with 
constant k inside a tube filled with viscous fluid with damping coefficient c. The motion 
of the particle is constrained parallel to the b3 axis only. At rest, the particle mass lies 
along the rotor’s bi coordinate axis, at a distance a from the rotor's center of mass. Three 
balancing masses, m2, m3, and m 4, each equal to the mass of the mass-spring-dashpot 
mechanism, are rigidly fixed a distance a on the b2, -bi, and -b2 axes. These masses 
maintain the axisymmetry of the system about the b3 axis. Displacement of the particle 
mass mi is denoted by the variable 2. A simplification in this paper of the Mingori dual- 
spinner system is the assumption of a massless spring-dashpot system. Thus the particle 
of the mass-spring-dashpot system, mi, is the same mass as the corresponding three other 
balancing masses, m2, m3, and m^. The platform can be described in a similar manner, 
with all notation modified with the prime mark. 

The dual-spin system center of mass coordinate axes rotates at the same rate as the 
rotor coordinate axes, but is located along the b3 axis at a distance z cm from the rotor center 
of mass. This distance will vary as the particle masses are displaced. Relating the dual- 
spin system's coordinate axes to the rotor coordinate axes was arbitrary. Equivalent results 
would be achieved by selecting the platform coordinate axes instead. 
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B. ANGULAR MOMENTUM 



The angular momentum of the Mingori dual-spin, quasi-rigid, axisymmetric system is 
now derived from first principles. 

The angular momentum (moment of momentum) of a particle of mass m; located in 
body B is defined as 

hi = r, x m, N \ l (108) 



where h, is the angular momentum of the i th particle with respect to the system center of 
mass, r, specifies the location of the i th particle with respect to the system center of mass. 



V_ dR i 



and V = 



dt 



N 



is the absolute velocity of the i th particle with respect to the inertial 



reference frame N. Expressing the absolute velocity in the system reference frame 



Nr»i N-cm N-B 

V = V +r,+ to xr,' 



(109) 



where B represents the reference coordinate axes of the system. For the following 
derivation, all displacements and velocities are referenced with respect to the rotor (b/) 
coordinate axes. The angular momentum vector is rewritten as 



. (N~cm . N~B \ 

h, = r, x m/ \ V + r,- + to x r,j 

Applying this equation to particle 1 with mass m\ and position 

H = a b! + 0 b 2 + (z — z cm ) b 3 



( 110 ) 



( 111 ) 
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one arrives at 



(112) 





( Z-ZcmY 0 -a(z-z cm ) 


1 ^ \ 




0 a 2 + (z - z cm f 0 


\ 0)2 l 




. -a(z-z cm ) o a 1 


\ O) 3 ) 





0 


-{z-z cm ) 


0 




' V cm x' 




0 






(Z - Zcm) 


0 


-a 




Vany 


| + 


~o{z — z cm ) 






0 


a 


0 




V cm z 




0 





For the rigid body, Equation (1 10) is applied to a differential volume at location 



r am = u bi + v b 2 + (w - z cm ) b 3 



Integrating, one arrives at 



h +M z} m 0 



h M = 



0 

0 



0 

— M z cm 

U o 



h +M z} m 
0 



M z cm 0 

0 0 

0 0 




(113) 



(114) 
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Similar calculations are performed for the remaining seven point masses and the platform. 
For the rotor, one arrives at 



h««4« = 

/ 



m 



h+M 2cm + 

(2a 2 + 4z c 2 « + 

' z 2 - 2 2 2cm 



- ma 2 



m 



h + M 2cm + 

(2 a 2 + 4 2cm +\ 
\ Z 2 - 2 2 2cm I 



-ma 2 
0 

73 + 4 ma 2 



0 ) l 
0)2 
0) 3 





0 






b, 


+ 


-ma 2 






b2 




. 0 






>3, 



0 (Af + 4m)(*») 
(M + 4 m)(z c «) 



-m2 



-m2 

0 



0 

0 



0 

0 

0 



V ) 

V cmy I 

Vcm. 



( 115 ) 
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Similarly, for the platform 



(116) 



h*r 

I 



♦ w = 

/,'+ M'{l-Z C mf + 

, f2 a' 2 + 4(1- z cm f + 
z' 2 + 2 z' (/ - z cm ) 



m 



m 



, (2 a' 2 + 4(1- z cm f + 
z' 2 +2 z'(l- Z C m) 



- m a z‘ cos (<r /) - m a’ z‘ sin (<r /) 



- m a.' z’ cos (o t) 



- m a,' z' sin (<r /) 



r / , i / /. 

h +4 m a 



0 )\ 
0)2 
0 ) 3 



\ 







— 1 






- m' a' <tcos(<t t) z’ + m' a' sin(a /) z' 






b ’l 


+ 


- m' a' a sin(<r /) z‘ - m' a ' cos(<7 1) z' 






b 2 




(h' + 4 m' a' 2 ) a 






bj / 



0 - (A/' + 4 m') (/ - Zcm) 

(A/' + 4 m')(l- z C m) 



- m 1 



0 

0 



0 

0 

0 






Vc, 



/ 



The angular momentum equations have terms corresponding to the velocity of the 
center of mass. Angular momentum, when taken about the center of mass of a system, by 
definition must be independent of the translation of the system's center of mass (with 
respect to an inertial reference frame). To verify that this occurs in the above equations, the 



equation for the position of the center of mass is substituted into the center of mass velocity 
terms, and then these equations are equated. If they have the same magnitude but opposite 
sign, they will cancel each other out and it will be proven that the system angular 
momentum is independent of the translation of the system center of mass. 
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One must prove 



-(■V. 1 (‘-4 JM 

Looking at the bi components, one finds 

- [(A/' + 4 m')(-l + z cm ) - m z] = [(M + 4 m) z cm - m z] 
Rewriting this, one obtains 

(m' + 4 m') l + m z + m z’ = {m + 4 m + M' + 4 m') z cm 
The center of mass with reference to the bi, b 2 , b 3 coordinate axes is 

V Mi z t + mi z -. , \ 

i m z + \M' + 4 m'j l + m' z' 



Zcm — 



(117) 



(118) 



(119) 



( 120 ) 



^ Mi + m- t M + 4 m +M' + 4 m 

i 

Substituting Equation (120) into Equation (1 19) results in 

{\t' + 4 m') / + mz + m' z' l [m * 4m m')j m ^ ' + ' Z ' l (121) 

\ M + 4 m +M' + 4 m J 



This reduces to 



(m' + 4 m') / + m z + m z = [m' + 4 m') / + m z + m z' (122) 

Therefore, the angular momentum along the bi axis is independent of the velocity of the 
center of mass. Verifying this along the other axes can be done in a similar manner. This 
proves that the angular momentum of the entire system about the system center of mass is 
independent of the translation of the system center of mass. The angular momentum of the 
entire system can be written as 

h = + 4 m + h W ' + 4 m ’ (123) 
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where 



hjK + 4«i = 



h + M Zen + 

(la 2 + 4 z} m +\ 

m \ z 2 -2z Zen I 



m 



-maz 



h + M Zen + 

fla 2 + 4 z} m +\ 

' Z 2 -lzZen I 
0 



0 

-maz 

0 



and 



/ 



♦ Am — 

h'+ M'{l-Zenf + 

,j2a' 2 + 4(l-Zenf + 

Z' 2 +2Z'(1-Zen) 

0 



m 



I 2 ’ + M'(l- Zcnf + 

, (2 a' 2 + 4 (l - Zenf + 
m ‘ Z' 2 +2Z’{1-Zen) 



-maz 

0 

/3 + 4 ma 2 



- m' a' z' cos (<T f) - m' a' z' sin (<r f) 



0)1 



a>2 1 



0)3 



( 124 ) 



- m' a' 2 ' cos (cr /) 



- m' a' z' sin (<r f) 



/ 3 ' + 4m a’ 2 



- m a’ o cos(<r t) z + m a' sin(<y t) 2 ' 
-m a' a sin(<r t) z' - m a’ cos(<7 t) z' 
(/ 3 ' + 4m' a' 2 ) <r 



( 125 ) 



0)1 



0)2 



0)3 
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C. ENERGY 



The kinetic energy, E, and the total energy, E to tai> are now derived. From first 
principles the kinetic energy of a differential particle is written as 

E<tm = \\ = £{( V") * ( N \ dm ) ) dm (126) 

where the absolute velocity is defined from before 



Nrrdm Nr'cm . tf-fl 

v = V + r dm + to xr dm 



(109) 



The kinetic energy of the Mingori system is determined by integrating the differential 
kinetic energy over the rotor and over the platform, and summing the differential kinetic 
energy equation over the eight point masses to arrive at 

4 

I 

i = 1 



£ = -L 
L 2 



J Rotor 1 — 1 

if !( v - ) ■ ( V' ) W* t i K V ) • ( V ) I m! 

Jn+n. r = 1 



Performing the steps on particle mj, the following is obtained 



(128) 



Ei = \rni J, 
^ i = l 



i(V m - V m Mh •r 1 )+((V m xri). (V m xri)) + 

2 (V" ■ r, ) + 2 (V" • ( N ST X ri) ) + 2 (r , ■ (V x n) ) 



Substituting in the values, one arrives at 



£ > = 2 m * 



v} m X + V2my+Vtmz + (z- Zcm? + (* “ Zcm? (£)\ + ’ 
d 2 0>3 + (z - Z cm ? (o\-2a(z- Zcm) CO\ CO3 + 
d CO2 + 2 (z V cm z — Z C m V cm z) 

2 / z ^2 ^cm j: — %cm ®2 V cm x + ^ ^3 V C m y ~ \ ^ 

Z CO\ V cm y + Z cm 0) 1 Vcm y ~ & ®2 V C m 2 
2(-azfi>2+0rcm<»2) 



(129) 
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This operation is then performed over the other three particles in the rotor. Performing 
these same steps on the body of the rotor, 

(130) 



Em = 




( N \ cm • N \ cm ) + (r 1 • n) + ( (^ cm x n ) . (N~cm x n)) + 

2 (V m • r ! ) + 2 (V"* • ( N a cm x r i) ) + 2 (f 1 • ( N n cm X r l) ) 



dM 



Substituting in values, and specifying the position vector of the differential particle as 

*dM = U bi + V b 2 + (w - z cm ) b 3 (131) 

then 



f 



Em = 



_1 



(132) 



V'cm x + Vcm y + V 'em z ^ Z cm + — 2 W Z cm CO 2 + 

z2 m 00% + V 2 CD* - 2 W V 0)2 CO 3 + 2 V Z cm 0)2 CO 3 + W 2 CO * + 
z} m col + u 2 col - 2 w z cm col - 2 u w co\ 0)3 + 2 u z cm <0i 0)3 + 



V 2 co\-2vu(0\ CO 2 + U 2 col + 2(-z cm V cmz ) + 

2 /vv CO2 V cmx~ Zcm CO2 ^cm x ~ v CO3 V cm x ~ w ® 1 V, 

( Z cm ©1 V cm y + U 0)3 V cm y + V <0i V c/n z - U fi> 2 1 
2 (- V Z cm <0] + U Z cm C0 2 ) 



dM 



cm y 

v, 



cm z 



Integrating over the limits of the symmetrical body, the kinetic energy of the rotor becomes 



Em = \ M 



V 2 

v cm x 
i 2 Z cm 



+ Vcm y + Vcm z + Z? m + z} m (<*> \ + 

I\ coj +I2CO2+ I3 CO3 - 
0)2 Vcm x + 2 Z cm C 0 \ V cm y ~ 2 V cm z j 



(133) 



These same steps are performed on the platform and on the remaining seven particles. 
These equations are then combined and simplified. 
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The total kinetic energy for the Mingori system is 



E = 



I 



(134) 

2 ~ 3 / total {y cm x + ^cm y + ^cm z) + ^ 

J (/i CD? + h CD? + /3 <0? ) + £(/i' 0)f + 12 col + V (0>3 + of) + 
l-(z 2 +z 2 {cd 2 + tofJ-Zaz o>i (©3) + 

^•(z' 2 +z' 2 (cd 2 + cd?) -2a' 1' cos(crr) coi C03 + 2 a' z' sin(ar) c 02 m3) + 

m(-flz m 2) + 

m'(-a' z' cos(ar) CD2 + a' z' sin(crr) cdi - a' z' crcos(crr) cdi - a' z' crsin(crr) CD2) 

I (cd 2 + CD 2 ) (m z + m' z') 2 - (m z + m z'f 

+ 2 Mp + 

/\2 






^-(cd 2 + cd 2 )(M + 4 m)/ 2 



M' + 4m' 

l mT 



^ (of + ^ 2 2 ) ( m ' + 4 m') / 2 ( MjMm j 2 ~ 

l z(cD 2 + CD?) |‘ 



mz + m' z' + (m' + 4 m') /) 



Mj 

, ,( 2 2 \im z + m' z' -(M + 4 m)l\ 

rz(co l+ co 2 )[ — ] 



The total energy of the Mingori system can be determined by adding the kinetic 
energy of the system to the potential energy of the system. The only potential energy of the 
system is the energy stored in the springs of the mass-spring-dashpot damper. From 
Hooke's law, the system potential energy is written as 

U (z) = l*z 2 +l*'z' 2 (135) 



The total energy of the Mingori system is then 

Etotai = E + U = £ + l*z 2 + ±*'z' 2 (136) 
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D. MINGORI’S EQUATIONS OF MOTION 

Mingori's equations of nx>tion [Ref 3] are written as follows 



A Q)\ -(A - C) ci>2 (&3 -J3 GCO2 + 2 Mt CCcoi + Mt C (©1 + 0>2 CO3) + 

( -2(C+/2)[z(^l-G>2 0 > 3 ) + ZG>l]+ l 

\z [-2 2 z co\ + z (fi)i - CO2 o)3)-a(d>3 + ce>i<a>2)] ) 

-2(C -li)[z'(<i)i- (02(03) + z' a>i]+ \ 

-2 £ (D\ + 2 z' (D\ + z'{(Di - 0)2(02) - a' cosy(q>2 + <*>16)2)] + / = C 
d f sin V [r + z'[{(03+af-(ol]) I 

A o>2 -(C - A) coi CO3 -J 3 gcoi + 2 M t CC^2 + Mt £ 2 (6)2 + a>\ 0)2) + 

I —2 (C / 2 ) [z (o>2 + 0)1 (t) 3 ) + z o> 2 ]+ I 

|z [-2 C Ct >2+ 2 z(0 2 + z {0)2 + <0\ 03)] -a[z + z (©3 - fi> 2 )] ) 

-2(C-/i)[z'fe + <0i<u 3 ) + z'a> 2 ] + 1 

m z'[-2 C 0)2 + 2 z 0)2 + z'(o)2 + 0)\0)3)]- = 0 

a' cosy/ {’z + z' [ {(02 + of -(of) )-a' smy/z'[d)2- o)\ 0)2) } 

C 6)2 - m a [2 z o)\ + z (©j - o) 2 (o 3 ) ] - 
rri a sin y /[2 z o)2 + z'(q ) 2 + q ) 1 0)2)]- 
m a cos y /[2 z 0 )i + z ( 6 ) 1 -o ) 2 a> 3 )] = 0 



m 



m (1 - p)z-m pz -ma[o>2- (0\ 0)2)- 
m ((of + (ofyz (1 -p)- 12 -p' z'] + c z + k z = 0 

-m p'z + m'( 1 - p') z' + 

m a' {sin y/[o)i + 0) 2 (0)3 + 2 a) ] - cos - G>i (0)3 + 2 <x)] } - 
m (o)f + (l -p') + li-pz] + c' z + k' z' = 0 



(137) 
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In the above equations, the following relationships are used 

Mj = M + M' + 4 mb + 4 mb 



AT + 4 mb_ 

Mj 




(138) 



Chapter IV discusses how these equations are adapted for use in the numerical integration 
routine. 

E. NUTATION ANGLE 

The nutation angle of the dual-spin system can be determined by substituting 
Equations (124) and (125) into Equation (34) to arrive at 



(139) 



6 - cos -1 





+ [li +4 m a' 2 ) <7 



IN 
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F. CORE ENERGY 



The platform core energy and the rotor core energy were defined in Equations (47) 
and (48). The nutation angle as a function of core energy was then developed and is stated 
in Equations (101) and (106). These equation can be used to predict the nutation angle of 
the system as a function of time. Because the core energy as a function of time is not 
available for the prediction, a postulated core energy must be developed. The initial value 
of the actual core energy and the postulated core energy must match and is dictated by the 
initial conditions. Parameters of the postulated core energy must be selected to accurately 
model the actual core energy, to include the final energy state and the rate at which the it 
approaches the final energy state. Two models were considered. 

1. Exponential Core Energy Model 

An exponential representation of rotor and platform core energy as a function of 
time are expressed as follows 

Ec postulated (?) = {ECo ~ Ec final) + Ec final (140) 

Ec postulated' M = {Ec 0 ' ~ E C final ) + Ec final' (141 ) 

The initial core energies, Ec 0 , Ec 0 '> are determined by the initial conditions of the system. 
The final core energies, Ec final* Ec final* as well as the exponential factors r, /, must be 
selected. The methodology for selecting these values is explained in Chapter V, Computer 
Analysis. 

2. Verhulst Logistic Core Energy Model 

The exponential model, as will be explained in Chapter V, has excellent 
agreement for stable conditions, but performs poorly for the unstable conditions. As an 
alternative model, the following first order differential model is selected for the rotor and 
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platform respectively. 



E C postulated (0 — 



E C postulated (0 — 



(gcJ i E C anal) 



Ec 0 + {Ec final - Ec 0 ) ^~ r ‘^ 

(EcfijEcfiJ) 



(142) 



(143) 



Ec :+(Ec final' -Eco')^ 1 ) 

This type of equation was first introduced by P. F. Verhulst to model human and 
other populations [Ref 7]. It is often referred to as the Verhulst equation or the logistic 
equation. Although population dynamics appears to be unrelated to the stability of a dual- 
spin system, the behavior over time of this equation compared to the stable and unstable 
systems provides some insight. Figure (4) shows the Verhulst equation versus time for 
varying initial conditions. 




44 



It can be noted that if the initial value of the variable (in our case core energy) is 
greater than the equilibrium value, it will approach the equilibrium value in a manner similar 
to that of an exponential decay. For initial values that are less than but within one half of 
the equilibrium value, the variable will asymptotically approach the equilibrium value. If 
the initial value is less than one half of the initial value, the variable over time has a 
somewhat different shape as it approaches equilibrium. Initially the slope is very small, but 
increases to a maximum at about the one-half of the equilibrium position. After this 
inflection point the variable approaches equilibrium asymptotically. The value of No equals 
K is an asymptotically stable equilibrium. A value of No equals zero is an unstable 
equilibrium. If the value of No is slightly greater than zero, then N(t) will, as t —»<*>, 
achieve the stable equilibrium of K. These types of curves will be utilized to describe both 
the stable and unstable dual-spin system. The initial core energies, Ec* Ec 0 '> would be 
determined by the initial conditions of the system. The final core energies, Eq final . 
Ec final* must be known or a best estimate used. The exponential factors, r, /, must be 
determined experimentally, and are a function of the system parameters and initial 
conditions. Chapter V provides additional details on the selection process. 
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IV. NUMERICAL SIMULATION 



A. NUMERICAL SIMULATION EQUATIONS OF MOTION 

The equations for the dynamical quantities needed for analysis of the dual-spin 
system were previously developed. Manipulation was required to make these equations 
suitable for numerical analysis, and is described below. The computer code is included as 
Appendix B. 

1. Mingori’s Equations of Motion 

The equations of motion for the dual-spin system are listed as Equation set ( 137 ) 
and ( 138 ). The Runge-Kutta numerical integration routine necessitated that the equations 
of motion be a set of first order differential equations. Mathematical manipulation was 
required to put them in this form. Variables representing groups of terms (A„ £„ C;, etc) 
were introduced to simplify the manipulations of the equations and provide a suitable 
format for incorporating the equations into the computer program. Mingori's equations can 
be expressed in terms of these variables and the five time-dependent variables of motion 

A26 &>i + A27 + A24 z' + A28 = 0 

&25 + ^23 C >3 + #13 Z + i?21 z' + #26 = 0 

C10 0 )\ + C5 Ct>2 + Cffl5 + Cu — 0 ( 144 ) 

D\ z + D2 z' + D3 o>2 + D% — 0 
£1 z + £2 z + Ei 6 >\ + £5 ci>2 + £10 = 0 

where the variables A t , Bi, C iy £>„ £„ F t , and Z, are defined in the Notation section. 
Clearly, these equations are highly coupled. Through a series of manipulations, including 
substitution and combining sets of equations to eliminate common variables, one can arrive 
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at a series of first order differential equations suitable for numerical integration techniques 



£ = * 
d t 



d/_ 

di 



= 2 



dz. 

dt 









Ml 
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7 Ql 7 . Ql 
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Z ‘ + H 
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( 145 ) 
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The equations remain coupled, so the sequence in which the equations are numerically 
integrated is important. Because the equations for 2 and the equation for 2' are functions of 
2 and 2', 2 and 2' must be integrated first. Similarly, 2 and 2' must be integrated before 0)3, 

Q) 3 before fi>i» and 6) 1 before 6>i. In the computer program, the seven variables that 
describe the motion are stored in a seven column matrix where fih = y[1][j], (O2 = y [2] [j], 

1 03 = y[3] 01, 2 = y [4] O'], 2 = y [5] 0], 2 ' = y[6] 01, and 2' = y(7] 01. The index j identifies the 
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matrix location for each set of saved variables over time. All other needed quantities can be 
calculated by using the seven time dependent variables describing the motion, and the 
parameters of the dual-spin system. These additional dynamical quantities are stored as one 
dimensional vectors in the computer. 

2. Angular Momentum 

Angular momentum was derived in Equations (123), (124), and (125). The 
angular momentum is a vector, but for verification of the conservation of angular 
momentum, only the magnitude is required. Therefore, the three vectorial components for 
both the rotor and the platform are computed individually as HI, H2, H3, Hip, H2p, H3p, 
and then the magnitude of the angular momentum, h[j], is determined and stored. 

3. Nutation Angle 

The nutation angle is determined using the previously determined relation 



6 - cos -1 



(h • b 3 

l N 



cos -1 



(h + h'\ 
{ N / 



cos 



-1 l^ s + I* 

N ; 



(34) 



Since it is a function of angular momentum, it can now be calculated and stored as 

theta [j]. 

4. Energy 

The total energy of the system, Equation (136), is determined by the sum of the 
kinetic energy, as derived in Equation (134), and the potential energy of the particle masses 
of the mass-spring-dashpot damper system. Equation (135). It is written as 



Etotai = E+U = £+l*z 2 + ±-*'z ' 2 ( 136 ) 

In the computer program the kinetic energy, ke [j], is the sum of the constituents of 
Equation (134), identified as T1 through T13. The potential energy is not explicitly stated 
in the program, but instead is included in the calculation of the total energy, etotai [j]. 
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5. Core Energy and Postulated Core Energy 

The platform and rotor core energy of the dual-spin system is defined in 
Equations (47) and (48) respectively. All of the variables are available, such that in the 
computer program the energies are readily calculated as Ecp [j] and Ec 01- The postulated 
nutation angle may now be computed. Equation (101) indicates a sign must be selected 
depending on whether the system is stable or unstable. The calculation of nutation angle 
over time in Equation (34) will indicate stability. The computer program incorporates 
conditional statements to assign the correct sign in the postulated nutation angle expression. 
Equations (101) and (106). In the derivation leading up to the postulated nutation angle, o 
has been defined as the relative rotation rate of the rotor with respect to the platform, and is 
a positive value. This provides a positive contribution to the angular momentum vector, 
thus providing stability to the system. In Mingori’s equations of motion, the reference 
coordinate axes are fixed to the rotor, resulting in the relative rotation of the platform in the 
counter-clockwise, or negative, direction. The postulated nutation angle equations were 
developed using the former reference frame. To compensate for the difference in the 
reference frames, the postulated nutation angle equations in the computer program have o 
replaced by -o. 

The postulated core energy as a function of time is computed for the platform and 
the rotor, using Equations (140) and (141) for the exponential model and Equations (142) 
and (143) for the Verhulst model. A postulated core energy model that accurately models 
the actual core energy required several iterations to find the proper values of the exponential 
factor. 

Once the postulated core energy models are computed, the Q' and Q values, as 
defined in Equations (100) and (107), are determined. The postulated nutation angles, 
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etah (j], and eta [j], are computed using Equations (101) and (106) with the proper sign 
previously determined. 

The revised stability criterion, Equation (97), requires the time rate of change of 
the core energy and the nutation frequency. Since both of these parameters may vary with 
time, an estimate is made to determine if the modified stability criteria correctly predicts 
stability or instability. The time rate of change of the core energy is approximated by taking 
the difference of the final and the initial core energy values, and dividing by the time of the 
simulation. This quantity is then divided by the initial nutation frequency. This stability 
quantity is computed for both the rotor and the core, but is indicated in Table (1), Summary 
of Analyzed Cases, simply as a negative or positive value. 

B. COMPUTER PROGRAM 

Essential to the verification of the stability criteria is the computer program that 
integrates the dual-spin system equations of motion, calculates other dynamical quantities 
of motion, and graphs the results. The system used was a Sun SPARC 2 workstation with 
the computer code written in C. Intermediate graphics results were created using an in- 
house computer graphics program. Final graphics output was performed by sending 
output data to the Deltagraph graphics package. The sequence of steps of the computer 
code are explained below. The computer code is included as Appendix B. 

1. Initialization 

The main computer program is compiled, along with the header file 'rkk.h' and the 
function 'derivs' immediately preceding it. The function 'derivs' contains Mingori's 
equations of motion rewritten into Equation (145), suitable for numerical integration. 
Included in the header file 'rkk.h' is the numerical integration routine ’rk4,' the adaptive step 
size function 'rkqc,' and the driver for the numerical integration routine 'odeint.' Also 
included are functions for creating and freeing the vectors and matrices used by the 
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computer to store the data, and a function for error messages. 

a. Variables 

The variables required for both the numerical integration routine and 
subsequent calculations are defined as global variables before the main program. Variables 
used only within a specific function or only in the main program are defined in their 
respective functions. Symbolic constants are also defined before the main program using 
the #define statement 

The time dependent variables of motion of the Mingori system, 
©i, © 2, 6)3, 2, 2, 2', 2', are stored in a seven column matrix. All other parameters to be 
graphed are stored in one dimensional vectors. Storage in the vectors and the matrix 
permits retrieval by the graphics subroutine for plotting the variable over time. 

b. Input File 

The main computer program scans the input file for the necessary parameters. 
All parameters of the Mingori system, l\, I2, h , M , m, a, k, c, l\\ I2, 1 3', a\ k\ 

c\ and L, are controlled with the input file. Also, the initial conditions for the time 
dependent variables of the system, ©1, ©2. £>3, z, z, z\ z', are specified. The length of 
time of the simulation and a variable determining the desired accuracy are specified. 
Finally, the exponential factor and the final energy for the postulated core energy functions 
are read. 

2. Preliminary Calculations 

After the input values are read in, initial calculations are performed prior to the 
numerical integration routine. All of the defined terms used by D. L. Mingori [Ref 3 ] in his 
non-linear equations of motion are computed, as well as definitions required for the core 
energy calculations. Is, Is > Is total * 1 1 , It > It total, <tnd h r igid. 

An option is available to specify critical damping in either the platform or the 
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rotor. By using the number 1000.0 in the input file for the damping coefficient, the main 
program will automatically compute the coefficient required for critical damping of the 
mass-spring-dashpot system and then use this value in all subsequent calculations. 

Also computed is the factor 'dxsav', used in determining when to save data. The 
steps between evaluating the equations of motion can become small, particularly when high 
accuracy is desired. The interval required for graphics resolution is not as restrictive. 
Accordingly, the variables are saved only if the step is greater than the previously saved 
step by the factor ’dxsav.’ 

3. Numerical Integration 

Mingori’s nonlinear equations of motion are numerically integrated by the fourth 
order Runge-Kutta method, with adaptive step size control. The computer code used is 
based on the Runge-Kutta method listed in Numerical Recipes in C, [Ref 8 ]. The adaptive 
step size control permits larger integration steps during smooth, well behaved portions of 
the functions, and smaller steps during the more irregular sections of the functions. The 
integration routine accuracy can be controlled by a variable in the input file. 

Mingori’s non-linear equations of motion are contained in the function 'derivs.' 
As described previously, the equations have been rewritten as a series of first order coupled 
differential equations suitable for numerical integration. 

The output of the time dependent variables of the system, co i, (Oi , co 3 , z, i, z\ z', 
is stored in an array of seven columns, with the number of rows required a function of the 
specified time interval and accuracy. A one dimensional vector is also created, with the 
time stored for each step saved. This permits plotting the time dependent variables and 
other quantities versus the time. 

4. Calculation of System Parameters 

With the array of the time dependent variables of motion as a function of time, the 
other system quantities listed in Section A may now be calculated, with the values stored in 
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vectors suitable for graphing. 

5. Graphics Output 

The remainder of the program is the necessary code for the in-house graphics 
program. The graphics program is initialized, and the graphics window is opened. The 
following parameters are then plotted over time: coi, co 2, 03-0)3 initial, 2, z, z', z', 

h ~ h initial , E ~ E initial, E-total ~ E- total initial , E c — E c initial, E c postulated ~ E c postulated initial, 
E c E c initial , Ec postulated Ec postulated initial • "H, tUld Tjf . The window is then closed 

and the initial conditions and pertinent time dependent variables and other quantities of the 
simulation are then printed. This displays the graphical representation of the behavior of 
the dual-spin system, and provides the actual initial and final values of the variables and 
associated quantities. 

For final graphics output, the Deltagraph graphics program is utilized. Computer 
simulation data is sent as an output file. The data is then manipulated into a suitable graphic 
with proper scaling and axes limits to best represent the dynamics of the computer 
simulation. The graphs are contained in Appendix A. 

6. Computer Program Validation 

The two aspects of the computer code requiring validation were the numerical 
integration routine and Mingori’s equations of motion with its associated dynamical 
quantities for the dual-spin system. 

The validation of the numerical integration routine was performed by using the 
equations of motion for a simple torque-free axisymmetric body. Various initial conditions 
were read in, and the dynamical quantities were then plotted versus time. The plotted 
behavior was then compared to the actual response of the system. 

The validation of Mingori’s equations of motion and the associated dynamical 
quantities required a more comprehensive approach. The equations of motion required 
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validation for proper behavior of each time-dependent variable of motion. Then the kinetic 
energy, total energy, and angular momentum must be verified. Several cases were run 
where the platform mass and inertia would become infinitesimally small. Then cases were 
run where the rotor mass and inertia would become infinitesimally small. In each of these 
cases, subcases were run where the mass-spring-dashpot system would be large, 
dominating the dynamics, to the subcase where it became very small, so the system 
approximates a rigid body. These different scenarios would uncouple and isolate the 
various parts of the equations of motion to verify proper derivation of the equations. In 
general, the platform and the rotor were tested under conditions varying from those in a 
rigid body scenario to those in a lightly damped body. The system was then tested as a 
rigid dual-spin system, as a dual-spin system with rotor damping only, and as a dual spin 
system with platform damping only. The dual-spin system was then tested with both the 
rotor and the platform containing dampers. For each case, the time dependent variables of 
motion were plotted and analyzed. In all cases, the angular momentum was compared with 
the initial angular momentum. Since all cases were torque free, the angular momentum 
must remain constant. In all cases explored, the angular momentum verified the 
correctness of the equations of motion. The code validation cases were not included in this 
thesis due to the large number of graphs and data that were required to establish validation. 
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V. ANALYSIS 



A. INTRODUCTION 

1. Objective 

The development of the revised energy-sink stability theory has been presented in 
Chapter II, and the equations of motion for a dual-spin system are contained in Chapter HI. 
With the numerical integration code of Chapter IV, the revised energy sink stability theory 
can now be verified. The core energy of the system is plotted as a function of time and 
compared to the total energy for agreement with theory. The stability criterion computed 
from the numerical simulation must then agree with Equation (97). Conservation of 
angular momentum is used in each case to verify the correctness of the equations, and to 
ensure the accuracy of the numerical integration routine. To approximate the core energies 
over time, an exponential model and a model based on logistic growth are then explored. 
The logistic growth model uses an equation presented by Verhulst [Ref 7], and is referred 
to as the Verhulst model. Postulated models of core energy and the associated nutation 
angles are then compared with the actual core energies and nutation angle for agreement 

2. Numerical Simulation Cases 

Four distinct cases needed to be addressed. With the inertia ratio of the dual-spin 
system greater than one, a stable case and an unstable case are analyzed (Cases (1) and 
(4)). With the inertia ratio less than one, a stable case and an unstable case are also 
analyzed (Cases (2) and (3)). For these four cases, an exponential model is used to 
postulate the core energy and the corresponding nutation angle. For the case of the inertia 
ratio greater than one and unstable, and the two cases of the inertia ratio less than one, 
stable and unstable, the Verhulst model is postulated (Cases (5), (6), and (7)). The case of 
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the inertia ratio greater than one and stable was not included; from the excellent agreement 
of Cases (5), (6), and (7), one can see that excellent agreement would also occur for this 
trivial case, making it unnecessary to include. All seven cases are summarized in Table (1). 
For each case, the dynamical quantities are plotted versus time until the system reaches a 
stable state. Each quantity is then plotted for the first one hundred seconds to show detail. 





h total 




ic 








CASE 


h total 


STABILITY 


X 


X' 


MODEL 


COMMENTS 


1 


1.039 


stable 


negative 


negative 


exponential 


large damping 
in platform 


2 


0.907 


stable 


negative 


negative 


exponential 


large damping 
in platform 


3 


0.897 


unstable 


positive 


positive 


exponential 


insufficient damping 
in platform 


4 


1.025 


unstable 


positive 


positive 


exponential 


insufficient damping 
in platform 


5 


0.907 


stable 


negative 


negative 


Verhulst 


same conditions 
as Case 2 


6 


0.897 


unstable 


positive 


positive 


Verhulst 


same conditions 
as Case 3 


7 


1.025 


unstable 


positive 


positive 


Verhulst 


same conditions 
as Case 4 



Table (1) Summary of Analyzed Cases 



In Appe~ "x A, tables (2) through (8) lis" ic system parameters, core energy 
parameters, and conditions. Immediately folk ing these tables are the graphs of the 
important dynar. ^ quantities. A typical geosynchronous dual-spin satellite was selected 
for the numerical simulation. The platform and rotor masses and inertias are listed as part 
of the initial conditions, and are the same for all cases. The inertia ratio is made greater 
than or less than one by selecting L, the distance between the rotor and platform centers of 
mass, to be 0.3 or 1.0 meter respectively. As a default set of values, all of the mass- 
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spring-dashpot systems have m and m equal to 1.0 kg, k and k' equal to 1.0 N/m, and c 
and c' equal to 1.0 kg/sec. To establish the stable cases, additional energy dissipation is 
required in the platform. To achieve this, cases (1), (2), and (5) have m increased to 20.0 
kg and c' increased to 10.0 kg/sec. In all seven cases, the initial conditions are the same. 
The platform rotates at a geosynchronous angular velocity, the rotor spins at the higher rate 
of 1.5 rad/sec, and a perturbation is introduced by an initial transverse angular velocity of 
0.1 rad/sec. The mass-spring-dashpot systems have no initial displacement or initial 
velocity. The core energy parameters are also listed in the tables. The initial core energies 
are determined by the system’s initial conditions. The final core energies were taken from 
the numerical simulation data. The exponential factors, r and /, were then determined 
through an iterative process to best fit the modeled core energy to the actual core energy. 

B. DISCUSSION 

The individual cases can now be analyzed. The graphs of the dynamical quantities 
are explored, and the data will affirm the revised stability theory. 

1. Angular Momentum 

For each case, the angular momentum is plotted versus time, Figures (5), (6), 
(16), (17), (27), (28), (38), (39), (49), (50), (60), (61), (71), and (72). Because the dual- 
spin system has no external forces, angular momentum must be conserved. For the stable 
cases, there is excellent agreement, with angular momentum varying by less than one one- 
hundredth of a percent over the length of the data run. This confirms the equations of 
motion and validates the accuracy of the numerical integration routine. For the unstable 
cases, the angular momentum percent difference increases to approximately four one- 
hundredths of a percent. This is attributed to the increased dynamics of the system as it 
establishes the spin about the transverse axis, introducing very small errors in the numerical 
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integration routine. The errors remain very small, and the equations of motion and the 
numerical integration routine re main accurate. 

2. Total Energy, Platform and Rotor Core Energy 

The graphs of energy are in Figures (7) through (10), (18) through (21), (29) 
through (32), (40) through (43), (51) through (54), (62) through (65), and (73) through 
(76). Several observations can be made. 

The total energy curve reaches equilibrium before both the platform and the rotor 
core energies do. For Case 2 and 5, the total energy appears to reach a maximum and then 
drop down to a final equilibrium value. Careful comparison of Figures (18) and (49), 
shows the same system with the same initial conditions, but with a slightly different curve. 
It can be deduced that the energy is not steady, but is still oscillating. The sampling 
frequency coincidentally saved values near the same magnitude of energy in the region 
from about 3000 seconds to 10000 seconds. 

The total energy for Cases (1), (3), (5), and (6) decreased. For Cases (2), (4), 
and (7), representing both stable and unstable systems, it increased. This can be explained 
easiest with the use of equations (26) and (27). For Cases (2), (4), and (7), the system is 
settling out about the axis with the minimum moment of inertia. Because there are no 
external forces, the angular momentum is constant In the equation for angular momentum. 
Equation (26), the inertia terms are squared. But in the equation for energy, Equation (27), 
the inertia terms are not squared. Therefore, as the angular velocity transfers to the axis 
with the minimum moment of inertia, Equations (26) and (27) show that the energy will 
increase. These equations apply to a simple dual-spin system. Equations (123) and (136) 
for the Mingori dual-spin system would show the same result, provided the energy 
absorbed by the mass-spring-dashpot system and the energy associated with the motor is 
less than the energy increase associated with transferring the spin to the axis of minimum 
moment of inertia. This is the situation for Cases (2), (4), and (7). 
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The total energy value will always be between that of the platform core energy 
and the rotor core energy. The platform core energy will have smaller values as its 
equation assumes that the rotor is rotating at the same rate as the platform, thereby not 
accounting for the large amount of kinetic energy associated with the rotor. The rotor core 
energy will be higher than the total energy as it assumes that the slowly spinning platform 
is rotating at the same rate as the rotor, providing the system with additional kinetic energy. 

The hundred-seconds graphs of total energy versus time shows curves with 
several different behaviors. The time rate of change of total energy includes the energy 
dissipation rate of the rotor and platform mass-spring-dashpot systems, and the rate of 
work due to the motor torque maintaining the relative rotation rate. In the hundred-seconds 
graphs, aside from the general slope of the curve, there is no apparent correlation between 
the specific behavior of the total energy to that of the core energies. Any relationship that 
exists is masked by the motor and mass-spring-dashpot system’s influences on total 
energy. 

3. Stability Criterion 

The revised stability criterion of equation (97) states that the time rate of change of 
the core energy over the respective nutation frequency must be less than or equal to zero. 
The sign of the stability criterion for each case must be determined. The time rate of change 
of the core energy was determined by dividing the final less the initial value of the core 
energy by the length of time of the case. The nutation frequencies were computed using 
initial conditions. The sign of the stability criterion is listed in Tables (1) through (8). The 
numerical value was not listed because the above assumptions used during the computation 
give it no merit. Cases with inertia ratios very close to one were intentionally selected to 
test the inertia ratio in this transition region. In all cases analyzed, the sign of the revised 
stability criterion was consistent with the stability of the system. 
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4. Postulated Core Energy and Nutation Angle 

An explicit relationship for core energy as a function of time does not exist If an 
equation describing core energy did exist, then by using Equations (101) and (106) the 
nutation angle as a function of core energy and time could be predicted for a dual-spin 
system. This could then be extended to an actual dual-spin satellite. Given a sufficient 
model of the satellite, the stability and the nutation angle as a function of time could be 
predicted. An objective of this thesis is to see if an equation for the core energy can be 
developed that would adequately describe the nutation angle as a function of time for both 
the stable and unstable conditions. The exponential model and the Verhulst logistic model 
were explored. 

a. Exponential Core Energy Model 

Observation of the core energy as a function of time for a stable system would 
lead one to conclude that it behaves in an exponential manner. Equations (140) and (141) 
are exponential representations of rotor and platform energy as a function of time. The 
initial core energy was determined by the initial conditions of the dual-spin system. The 
final core energy was determined by running the numerical simulation and calculating it at 
the end of the simulation. If this was not available, a value could be estimated. Applying 
the principle of conservation of angular momentum and noting that the system will 
eventually spin about one of the primary axes, then the equations for angular momentum 
and total energy can be used to solve for the angular velocity about that axis. Substituting 
this into Equation (47) or (48), one would arrive at an estimated final core energy. 
Although it would not be the actual final core energy, as motor torque contribution and the 
mass-spring-dashpot system’s energy dissipation was not accounted for, it would be 
sufficiently close to satisfy the computational requirements. The exponential factor is 
dependent upon the system parameters and initial conditions. For the cases presented, 
different values were tried until good agreement was established with the core energy 
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curve. Additional research would be required to determine a suitable exponential factor for 
an actual satellite, taking into account the parameters and the initial conditions. Cases 1 
through 4 contain the numerical simulation data for the exponential model. The actual and 
postulated curves for core energy are contained in Figures (11) through (15), (22) through 
(24), (33) through (37), and (44) through (48). For cases 1 and 2, both stable, there is 
excellent agreement between the core energy and the postulated core energy. This in turn, 
results in excellent agreement between the actual nutation angle and the modeled nutation 
angle. The actual nutation angle is determined using the angular momentum quantities, as 
shown . Equation (34). The time dependent variable required to compute the nutation 
angle is the angular velocity about the spin axis. The modeled nutation angle is determined 
in Equations (101) and (106), and is a function of only one time-dependent variable, core 
energy. Herein lies the potential of the core energy theory. A sufficient model of core 
energy over time, as in Cases 1 and 2, will provide an excellent prediction of nutation 
angle, without requiring any knowledge of the specific angular velocities of the system as a 
function of time. 

Cases 3 and 4 illustrate the exponential model for an unstable dual-spin 
system. The exponential model for core energy, and its associated nutation angle, rapidly 
approach their final values. The actual core energy and nutation angle, however, behave 
quite differently. The initial conditions have the system near an unstable equilibrium. The 
system moves from the unstable to the stable equilibrium slowly at first. It then increases 
the rate at which it approaches stable equilibrium, passes through an inflection point, and 
then approaches equilibrium asymptotically. It is clear that the exponential model 
represents this behavior poorly. Verhulst’s logistic equation was then addressed to 
determine its adequacy in modeling the dual spin system. 
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b. Verhulst Logistic Core Energy Model 



The Veihulst logistic equation was introduced in Chapter HI. Figure (4) and 
the associated description explains the characteristics of the equation. The general shape of 
the curves in Figure (4) is very similar to the stable and unstable cases of the dual-spin 
system. Case 5 is the Verhulst model of the Case 2 stable system. Figures (55) through 
(59) show that the Verhulst logistic equation can achieve excellent agreement with core 
energy and nutation angle, just as the exponential model did. 

Cases 6 and 7 are the same as Cases 3 and 4, except the Verhulst logistic 
equation is used to model the core energies. Figures (66) through (70) and (77) through 
(81) illustrate the modeled and actual core energies. Although the Verhulst model for rotor 
core energy had excellent agreement, it performs poorly when modeling the unstable 
system. The rotor core energy begins at an initial value and then decreases to an 
equilibrium value. Referencing Figure (4), the Verhulst logistic equation will model the 
rotor core energy exponentially. 

To take advantage of the Verhulst logistic equation’s curve beginning near the 
the unstable equilibrium and its progression to the asymptotically stable equilibrium, the 
core energy must increase over time. The platform core energy behaves in this manner as 
the nutation angle goes from a small angle to ninety degrees. Figures (68) through (70) 
and (79) through (81) show the Verhulst modeled and the actual platform core energies, 
and the modeled and the actual nutation angles. The agreement between the actual and 
modeled energies and nutation angles was very good. The general shape of curve was 
consistent, with only a slight deviation in the center of the curve, and then a small deviation 
as the nutation angle approaches the equilibrium value of ninety degrees. With this 
agreement, it is established that a core energy model exists that can represent both stable 
and unstable cases. For each of the cases, core energy is plotted versus nutation the angle, 
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Figures (15), (26), (37), (48), (59), (70), and (81). The graphics package permitted only a 
semi-log plot of rotor core energy. The plot did not provide any insight on relationships 
for the dynamics of the dual-spin system. For both the stable and unstable cases, the log- 
log plot between platform core energy and the nutation angle is linear from ninety degrees 
until about two degrees. For nutation angles less than two degrees, the platform core 
energy asymptotically approaches the equilibrium value. Case 5 and 6 is the same system 
with the same initial conditions, with the exception of the increased damping system in the 
platform for Case 5. The initial platform core energies and nutation angles are very nearly 
the same for each case, as shown in Figures (57), (58), (59), (68), (69), and (70). One 
can conclude by this observation and the definition of nutation angle as a function of core 
energy. Equation (101), that there exists a continuous platform core energy versus nutation 
angle curve. The appearance would look similar to the curve created by splicing Figures 
(59) and (70) together. By varying one parameter, for example platform energy dissipation 
rate, the system would progress along this curve and achieve equilibrium with a zero 
degree nutation angle or achieve equilibrium with a ninety degree nutation angle. This is 
what was done with cases 5 and 6. Adding as the third dimension to the curve the time rate 
of change of core energy would then reveal the stability, the initial direction, and the rate at 
which the system will arrive at the equilibrium condition. 

C. FURTHER RESEARCH 

The revised stability criterion for a dual-spin, quasi-rigid, axisymmetric system was 
established. Numerical simulation was then used to verify the revised stability theory. 
Further research could be conducted in several areas. The specific contributions to the total 
energy could provide some insight How much energy and in what manner does the motor 
torque contribute to the total energy for both the stable and the unstable cases? Also, by 
plotting the energy dissipation system’s contributions over time, one could determine its 
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effect on total energy, as well as comparing it to the core energy theory in Equation (97). 

Relationships for core energy as a function of time were explored. It was established 
that the Verhulst logistic equation could be applied to a stable dual-spin system. By 
changing its parameters, this equation also applied to the unstable dual-spin system with 
very good agreement. Further research could be directed to fmd one equation that would 
address both the stable and unstable cases of the dual-spin system without changing its 
parameters. The equation for logistic growth with a threshold [Ref. 7] as applied to the 
platform core energy shows promise in this regard. 



Given the initial conditions, the platform core energy will typically start at some 
intermediate value and will either increase or decrease as the dual-spin system reaches 
equilibrium with a nutation angle of either zero degrees or ninety degrees. Equation (146) 
would be well suited to model the platform core energy. There exists a platform core 
energy value, Eq unstable > where the system is at an unstable equilibrium. Equation (146) 
shows that at exactly this value, the rate of change of the platform core energy will be zero. 
For any value below the unstable equilibrium value, the core energy will approach the value 
of zero. For Case 5, the final core energy for the stable condition was 0.122 J. Although 
this final condition cannot be represented in Equation (146), it is sufficiendy close that the 
equation may still be used to represent the platform core energy. For any value above the 
unstable equilibrium value, the core energy will approach the equilibrium value of Eq final* 
associated with the nutation angle at ninety degrees. Initial conditions would determine the 
initial core energy. The final core energy can be estimated as described earlier in this 
chapter. Finally, the exponential factor / may then be determined based on the system 
parameters and initial conditions. 

Additional analysis can be performed on the revised stability criterion. The time rate 




(146) 
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of change of core energy over nutadonal frequency may be calculated and plotted versus 
time or versus other parameters. The behavior of this stability criterion and the magnitude 
of it for various conditions could provide some insight on postulating a relationship for 
core energy as a function of time. 
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VL CONCLUSION 



The existing energy-sink stability criterion was introduced as 




( 41 ) 



X X' 



It was shown that an inconsistency in the development disproves the existing 
assumption that the motor energy input exacdy balances shaft frictional losses. A revised 
energy-sink stability criterion was then developed based on Hubert’s definition of core 
energy and was presented as 



This criterion compliments the existing theory. Numerical simulation was required to 
validate the theory. The Mingori dual-spin, quasi-rigid, axisymmetric system was selected 
for the numerical simulation. Several cases were analyzed to verify the revised energy-sink 
stability criterion. By correctly postulating the platform or the rotor core energy, the 
stability of the system could be determined. Specific knowledge of the energy dissipation 
rates for both the platform and the rotor are no longer required. 

An exponential model and the Verhulst logistic model for core energy, and their 
relationship to the nutation angle, were explored. The exponential model had excellent 
agreement with the stable cases, but was inadequate in representing the unstable cases. The 
Verhulst logistic model established that an explicit relationship for core energy could be 
developed. Nutation angle as a function of core energy for the dual-spin system could then 
be predicted. The excellent agreement of the postulated core energy and nutation angle with 
the actual core energy and nutation angle confirms the revised energy-sink stability 
criterion. Additional research is required to find an optimum equation for core energy as a 
function of time to represent both the stable and unstable cases. 




(97) 



X X' XX' 
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APPENDIX A NUMERICAL SIMULATION DATA 



CASE 1: y- > 1, Stable, Exponential Model 



System Parameters 
Platform Rotor 



/]'= 1600 kgm 2 


/ 1 = 


1000 kgm 2 


1 2 = 1500 kgm 2 


/ 3 = 


1200 kgm 2 


M’ = 1000.0 kg 


M = 


700.0 kg 


m = 20.0 kg 


m = 


1.0 kg 


a = 1.0 m 


a = 


1.0 m 


k'= 1.0 

m 


k = 


i.o£ 

m 


c' = 10.0 |§- 
sec 


c = 


1.0 

sec 


L = 0.3 m 


Is total 
h total 


= 1.039 



Initial Conditions 
Platform Rotor 



2 ' = 


0.0 m 


2 = 


0.0 m 


z = 


0.0 ^ 
s 


2 = 


0.0 Ip- 
s 


003 ' = 


7.27 x l0- 5Qd - 
sec 


0)2 = 


1.5Qi 

sec 


0)1 = 


o.io^i 

sec 


0)2 = 


0.0 ^ 
sec 



Core Energy Parameters 
Platform Rotor 



E c initial ~ 


13.402 J 


E c initial “ 


3145.4 J 


Ec final ~ 


0.072 J 


E c final — 


3161.7 J 


II 

-U5* | ^ 


negative 


n 


negative 


/ = 


-.00134 s- 1 


r = 


-.00134 s- 1 



Table (2) Case 1 Parameters 
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Figure (5) Total Energy and Percent Difference Angular Momentum Versus Time 
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Figure (7) Total Energy and Rotor Core Energy Versus Time 



Case 1: Is/It>l, Stable, Exponential Model 
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Figure (8) Total Energy and Rotor Core Energy First 100 Seconds 
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Figure (9) Total Energy and Platform Core Energy Versus Time 




Figure (10) Total Energy and Platform Core Energy - First 100 Seconds 
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Case l:Is/It>l, Stable, Exponential Model 
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Figure (11) Modeled and Actual Rotor Core Energy and Nutation Angle Versus Time 



Case 1: Is/It>l, Stable, Exponential Model 
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Figure (12) Modeled and Actual Rotor Core Energy and Nutation Angle- First 100 Seconds 
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Case 1: Is/It>l, Stable, Exponential Model 
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Figure (13) Modeled and Actual Platform Core Energy and Nutation Angle Versus Time 
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Figure (15) Core Energy Versus Nutation Angle 
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CASE 2: ^ < 1 

It 



, Stable, Exponential Model 



System Parameters 
Platform Rotor 



(\ = 1600 kgm 2 


/ 1 = 


1000 kgm 2 


h'= 1500 kgm 2 


/ 3 = 


1200 kgm 2 


tf'= 1000.0 kg 


A/ = 


700.0 kg 


7t = 20.0 kg 


m = 


1.0 kg 


a = 1.0 m 


a = 


1.0 m 


k'= 1.0 £ 

m 


k = 


1.0 g. 

m 


c — 10.0 

sec 


c = 


1.0 M 

sec 


L= 1.0 m 


li total 

It total 


= 0.907 



Initial Condiut os 





Platform 


Rotor 




z' 


= 0.0 m 


z = 


O 

O 


m 


z' 


= 0.0 “1 

s 


z = 


O 

d 


m 

s 


«3' 


= 7.27 x lO" 5 

sec 


0 >. 3 = 


1.5 


rad 

sec 


©1 


= o.io 

sec 


6)2 = 


0.0 


rad 

sec 



Core Energy Parameters 
Platform Rotor 



E c initial “ 


15.34 J 


initial ~~ 


3147.3 J 


E c final = 


0.123 J 


E c final = 


3170.9 J 


II 


negative 


HS 5, 

II 


negative 


/ = 


-.00102 s- 1 


r = 


-.00102 s- 1 



Table (3) Case 2 Parameters 
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Figure ( 1 6) Total Energy and Percent Difference Angular Momentum V ersus Time 



Case 2: Is/It<l, Stable, Exponential Model 
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Figure (17) Total Energy and Percent Difference Angular Momentum - First 100 Seconds 
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Figure ( 1 8) Total Energy and Rotor Core Energy Versus Time 
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Figure (20) Total Energy and Platform Core Energy Versus Time 
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Figure (22) Modeled and Actual Rotor Core Energy and Nutation Angle Versus Time 
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Figure (23) Modeled, Actual Rotor Core Energy and Nutation Angle - First 100 Seconds 
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Figure (24) Modeled and Actual Platform Core Energy and Nutation Angle Versus Time 



Case 2: Is/It<l, Stable, Exponential Model 
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Figure (25) Modeled, Actual Platform Core Energy and Nutation Angle- First 100 Seconds 
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Figure (26) Core Energy Versus Nutation Angle 
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CASE 3: ^ < 1 

h 



, Unstable, Exponential Model 



System Parameters 
Platform Rotor 



/ t'= 


1600 kg m 2 


/ 1 = 


1000 kgm 2 


h' = 


1500 kgm 2 


/ 3 = 


1200 kgm 2 


M' = 


1000.0 kg 


M = 


700.0 kg 


m = 


1.0 kg 


m = 


1.0 kg 


a = 


1.0 m 


a- 


1.0 m 


k' = 


1.0 S 

m 


k = 


1.0 & 
m 


c' = 


1.0 — 
sec 


c = 


1.0 M. 

sec 


L - 


1.0 m 


I s total 
h total 


= 0.897 



Initial Conditions 
Platform Rotor 



i - 0.0 m 


z = 


0.0 m 


z' = 0.0 m 

s 


z = 


00 f 


o> 3 ' = 7.27xl0- 5 g| 


0>3 = 


i s rad. 

° sec 


o)i = 0.10 gd- 

sec 


tt>2 = 


0.0 sd 

sec 



Core Energy Parameters 
Platform Rotor 



E c initial — 


15.1 J 


Ec initial ~ 


3061.6 J 


E c final ~ 


1 160.9 J 


Ec final ~ 


1489.9 J 


* h"! 

II 


positive 


Ec 

X 


positive 


r' = 


-.0005 s -1 


r = 


-.0005 S" 1 



Table (4) Case 3 Parameters 
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Figure (27) Total Energy and Percent Difference Angular Momentum Versus Time 




Figure (28) Total Energy and Percent Difference Angular Momentum - First 100 Seconds 
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Figure (29) Total Energy and Rotor Core Energy Versus Time 
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Figure (30) Total Energy and Rotor Core Energy - First 100 Seconds 




Figure (3 1 ) Total Energy and Platform Core Energy Versus Time 
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Case 3: Is/It<l, Unstable, Exponential Model 
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Figure (32) Total Energy and Platform Core Energy - First 100 Seconds 
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Figure (33) Modeled and Actual Rotor Core Energy and Nutation Angle Versus Time 
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Figure (34) Modeled, Actual Rotor Core Energy and Nutation Angle - First 100 Seconds 
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Figure (35) Modeled and Actual Platform Core Energy and Nutation Angle Versus Time 
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Figure (36) Modeled, Actual Platform Core Energy and Nutation Angle- First 100 Seconds 
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CASE 4: ^ > 1 

It 



, Unstable, Exponential Model 



System Parameters 
Platform Rotor 



/l'= 


1600 kgm 2 


/ 1 = 


1000 kgm 2 


h' = 


1500 kgm 2 


/3 = 


1200 kgm 2 


M' - 


1000.0 kg 


M = 


700.0 kg 


m = 


1.0 kg 


m = 


1.0 kg 


a' = 


1.0 m 


a = 


1.0 m 


k' = 


i.o £ 

m 


k = 


i.o £ 

m 


c' = 


1.0 ^ 
sec 


c = 


1.0 

sec 


L - 


0.3 m 


I s total 
1 1 total 


= 1.025 



Initial Conditions 
Platform Rotor 



z - 0.0 m 


2 = 


0.0 m 


z = 0.0 


2 = 


00 f 


© 3 ' = 7.27 xio- 5 ^- 


a>3 = 


15 s 


©i = 0.10 ^ 

sec 


Ct>2 = 


0.0^ 



Core Energy Parameters 
Platform Rotor 



initial ~ 


13.2 J 


fc initial ~ 


3059.7 J 


E c final ~ 


1247.1 J 


E c final = 


1552.5 J 


Ec_ 

X 


positive 


II 

•&l 


positive 


/ = 


-.0005 s- 1 


r = 


-.0005 s- 1 



Table (5) Case 4 Parameters 
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Figure (38) Total Energy and Percent Difference Angular Momentum Versus Time 



Case 4: Is/It>l, Unstable, Exponential Model 
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Figure (39) Total Energy and Percent Difference Angular Momentum - First 100 Seconds 
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Figure (4 1 ) Total Energy and Rotor Core Energy - First 100 Seconds 
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Figure (42) Total Energy and Platform Core Energy Versus Time 
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Case 4: Is/It>l, Unstable, Exponential Model 
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Figure (43) Total Energy and Platform Core Energy - First 100 Seconds 
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Figure (46) Modeled and Actual Platform Core Energy and Nutation Angle Versus Time 



Case 4: Is/It>l, Unstable, Exponential Model 
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Figure (47) Modeled, Actual Platform Core Energy and Nutation Angle- First 100 Seconds 
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CASE 5: ^ < 1 

It 



, Stable, Verhulst Model 



System Parameters 



Platform 




Rotor 


I\ = 1600 kgm 2 


/i = 


1000 kgm 2 


I 3 '= 1500 kg m 2 


/3 = 


1200 kgm 2 


M' - 1000.0 kg 


M = 


700.0 kg 


m = 20.0 kg 


m = 


1.0 kg 


d = 1.0 m 


a = 


1.0 m 


k'= 1.0 

m 


* = 


i.o a 

m 


c — 10.0 

sec 


c = 


1.0 !§: 
sec 


L = 1.0 m 


Is total 
I 1 total 


= 0.907 



Initial Conditions 
Platform Rotor 



z' = 0.0 m 


Z = 


0.0 m 


z' = 0.0 n 


z - 


0.0 m 


<o 3 ’= 7.27x10-5^ 


(d = 


15 Ed. 

sec 


0)1 = 0.10 gl 


(O 2 = 


0.0 gd. 
sec 



Core Energy Parameters 
Platform Rotor 



E c initial ~ 


15.341 J 


Ec initial — 


3147.3 J 


Eq final ~ 


0.122 J 


E c final = 


3170.9 J 


II 


negative 


hs* 

II 


negative 


r' = 


-.0002 s- 1 


r = 


-.0002 s-> 



Table (6) Case 5 Parameters 
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Figure (49) Total Energy and Percent Difference Angular Momentum Versus Time 




Figure (50) Total Energy and Percent Difference Angular Momentum - First 100 Seconds 
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Case 5: Is/It<l, Stable, Verhulst Model 
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Figure (51) Total Energy and Rotor Core Energy Versus Time 



Case 5: Is/It<l, Stable, Verhulst Model 
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Figure (52) Total Energy and Rotor Core Energy - First 100 Seconds 
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Figure (55) Modeled and Actual Rotor Core Energy and Nutation Angle Versus Time 
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Case 5: Is/It<l, Stable, Verhulst Model 




Figure (57) Modeled and Actual Platform Core Energy and Nutation Angle Versus Time 



Case 5: Is/It<l, Stable, Verhulst Model 




Figure (58) Modeled, Actual Platform Core Energy and Nutation Angle- First 100 Seconds 
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Case 5: Is/It<l, Stable, Verhulst Model 




1 

nutation angle (degrees) 



100 



10 



>> 

8 ? 

<D 

c 

% 

o 

u 



0.1 



0.1 



Figure (59) Core Energy Versus Nutation Angle 
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CASE 6: ^ < 1 

It 



, Unstable, Verhulst Model 



System Parameters 
Platform Rotor 



Initial Conditions 
Platform Rotor 



r i'= 1600 kgm 2 


/ 1 = 


1000 kgm 2 


r 3 '= 1500 kgm 2 


/3 = 


1200 kgm 2 


r = 1000.0 kg 


M = 


700.0 kg 


n = 1.0 kg 


m = 


1.0 kg 


a = 1.0 m 


a = 


1.0 m 


t'= 1.0 g 

m 


* = 


1.0 g 

m 


:' = 1.0 M 

sec 


c- 


1.0 

sec 


L= 1.0 m 


li to'ai 

h total 


= 0.897 



z = 0.0 m 


z - 


0.0 m 


z' = 0.0 “■ 

s 


z = 


0.0 m 
s 


0)3 = 7.27 x 10- 5 

sec 


c 03 = 


1.5 J2d 

sec 


coi = 0.10 

sec 


0)2 = 


0.0 ^d 

sec 



Core Energy Parameters 
Platform Rotor 



E c initial ~ 


15.1 J 


E c initial — 


3061.6 J 


Ec final ~ 


1161.0 J 


E c fmal = 


1489.9 J 


Ec_ 

X 


positve 


II 


positive 


r' = 


-1 x 10- 6 s- 1 


r = 


-lx 10 -6 s -1 



Table (7) Case 6 Parameters 
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Case 6: Is/It<l, Unstable, Verhulst Model 
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Figure (60) Total Energy and Percent Difference Angular Momentum Versus Time 




Figure (61 ) Total Energy and Percent Difference Angular Momentum - First 100 Seconds 
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Figure (62) Total Energy and Rotor Core Energy Versus Time 
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Figure (63) Total Energy and Rotor Core Energy - First 100 Seconds 




Figure (64) Total Energy and Platform Core Energy Versus Time 
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Case 6: Is/It<l, Unstable, Verhulst Model 
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Figure (65) Total Energy and Platform Core Energy - First 100 Seconds 
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Case 6: Is/It<l, Unstable, Verhulst Model 
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Figure (66) Modeled and Actual Rotor Core Energy and Nutation Angle Versus Time 




Figure (67) Modeled, Actual Rotor Core Energy and Nutation Angle - First 100 Seconds 
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Case 6: Is/It<l, Unstable, Verhulst Model 
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Figure (68) Modeled and Actual Platform Core Energy and Nutation Angle Versus Time 



Case 6: Is/It<l, Unstable, Verhulst Model 




Figure (69) Modeled, Actual Platform Core Energy and Nutation Angle- First 100 Seconds 
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Figure (70) Core Energy Versus Nutation Angle 
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CASE 7: ^->1 

It 



, Unstable, Verhulst Model 



System Parameters 
Platform Rotor 



Initial Conditions 
Platform Rotor 



/l'= 


1600 kg m2 


/ 1 = 


1000 kgm 2 


h' = 


1500 kgm 2 


h = 


1200 kgm 2 


M' = 


1000.0 kg 


M = 


700.0 kg 


m = 


10 kg 


m = 


1.0 kg 


a' - 


1.0 m 


a = 


1.0 m 


k' = 


1.0 * 
m 


* = 


1.0 & 


c' = 


1.0 is- 

sec 


c = 


1.0 is. 

sec 


L = 


0.3 m 


Is total 


= 1.025 






1 1 total 





z = 0.0 m 


z = 


0.0 m 


z' = 0.0 ^ 


z = 


00“ 


(03 = 7.27 xlO- 5 ^ 


G^ = 


15 si 
sec 


(Oi = 0.10 Sd 

* sec 


0)2 = 


0.0 ^d 
sec 



Core Energy Parameters 
Platform Rotor 



Ec initial ~ 


13.2 J 


E c initial ~ 


3059.7 ; 


E c final = 


1247.1 J 


Ec final = 


1552.5; 


II 


positve 


II 


positive 


r' = 


-l.lxlO^s- 1 


r = 


-1.1 xlO- 6 



Table (8) Case 7 Parameters 
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Figure (71 ) Total Energy and Percent Difference Angular Momentum Versus Time 
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Figure (72) Total Energy and Percent Difference Angular Momentum - First 100 Seconds 
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Case 7: Is/It>l, Unstable, Verhulst Model 
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Figure (74) Total Energy and Rotor Core Energy - First 100 Seconds 
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Figure (75) Total Energy and Platform Core Energy Versus Time 




Figure (76) Total Energy and Platform Core Energy - First 100 Seconds 
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Case 7: Is/It>l, Unstable, Verhulst Model 
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Figure (77) Modeled and Actual Rotor Core Energy and Nutation Angle Versus Time 




Figure (78) Modeled, Actual Rotor Core Energy and Nutation Angle - First 100 Seconds 
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Figure (79) Modeled and Actual Platform Core Energy and Nutation Angle Versus Time 
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Case 7: Is/It>l , Unstable, Verhulst Model 
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Figure (80) Modeled, Actual Platform Core Energy and Nutation Angle- First 100 Seconds 
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APPENDIX B COMPUTER PROGRAM CODI 



/*****************•«•** declarations ***************** * 444 4 * / 

♦include <math.h> 

♦include <malloc.h> 

♦include <stdio.h> 

♦define MAXSTP 150000 
♦define TINY 1.0e-40 

♦ define PGROW -0.20 /* krqc */ 

♦define PSHRNK -0.25 /* krqc */ 

♦define FCOR 0.06665666666666666666666666666667 /* 1.0/15.0 krqc 4 / 
♦define SAFETY 0.9 /* krqc 4 / 

♦define ERRCON 6.0e-4 /* krqc •/ 

float **y-0, *xx-0; /* defining declaration rkdumb 4 / 

int kmax-0, kount-C; (* defining declaration odeint */ 

float *xp-0, * 4 yp-0, dxsav-0; /* defining declaration odeint */ 



/ 4 ode in? 4 / 
/ 4 ode i nt 4 / 



/* **************** * error function **********************/ 
void nrerror (error_text) 
char error_text [ ] ; 

{ 

void exit ( ) ; 

fprintf (stderr, "Numerical Recipes run-time error ... \n") ; 
fprintf (stderr, " *s\r." , error_text ) ; 
fprintf (stderr , " . . . now exiting to system. .. \n") ; 
exit (1) ; 

) 

/**************•***•»** vector function ************************/ 
float * vectorr (nl, nh/ 
int nl, nh; 

{ 

float *v; 

v - (float *) malloc ( (unsigned) (nh-nl + 1) *sizeof (float) ) ; 
if ( ! v) nrerror ( "allocation failure in vector ()") ; 
return v-nl; 

) 

/•************«••**«• matrix function *************************** 
float * *matrix (nrl , r.rh, nci, nch) 
int nrl, nrh, ncl, nch; 

{ 

int i; 
float * *m; 



it - (float **) malice { (unsigned) (nrh-nr 1 + 1 ) *sizeof (f loat 4 ) ) ; 
if (!m) nrerror ( "aliccation failure 1 in matrix ()"); 
m — nrl; 

for (i-nr 1; i<-nrh; i-**-) 

( 

m[i) - (float *) malloc ( (unsigned) (nch-ncl+1 ) * sizeof ( f loot ) ) ; 
if (!m(i]J nrerror ( "allocation failure 2 in matrixO"); 
m(i) — r.cl; 

) 

return m; 

J 

/*****•«•*•*«****•**• « ree vector function ********************/ 
void f ree^vector (v f nl , nh) 
float *v;’ 
int nl, nh; 
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free ( (char* ) (v + nl) ) ; 

) 

/**************♦********* free matrix function ************* 
void f ree_matrix (m, nrl, nrh, ncl, nch) 
float **m; 

int nrl, nrh, ncl, nch; 

{ 

int i: 

for (i-nrh; i>-nrl; i — ) free ((char*) (m[i)+ncl) ) ; 
free ((char*) (m+nrl) ) ; 

) 

/******♦***♦*********** r k 4 function ********************** / 
void rk4 ty, dydx, n,x,h,y out, derivs) 

float yl)/ dydx(], x, h, yout (] ; 
void ( *derivs ) ( ) ; 
int n; 

( 

int i; 

float xh, hh, h6, *dym, *dyn, *dyt, *yt, *vectorr(); 
void free vector () ; 



dym 


- vectorr (1 


, n) 


dyn 


- vectorr (1 


, n) 


dyt 


- vectorr(l 


, n) 


yt - 


vectorr (1, 


n) ; 


hh - 


h * 0 . 5 ; 




h6 - 


h/6.0; 




xh - 


x + hh ; 





for ( i— 1 ; i<-n; i++) yt[i) - y ( i ) +hh*dydx [ i ) ; /* first step •/ 

(•derivs) (xh, yt, dyt , dydx) ; 

for (i-1 ; i<-n; i++ ) yt(i] - y[i)+hh*dyt[i); 

(*derivs) (xh, yt , dym, dyt ) ; 
for (i-1 ; i<-n; i++ ) 

( 

yt I i ) - y ( i ) +h*dym ( i ) ; 
dyn[i) - dyt [ i) +dym (i) ; 

) 

(*derivs) (x + h, yt, dyt, dym) ; 

for (i-1 ; i<-n; i++ ) yout[i] - y ( i ] +h6* (dydx ( i ] +dyt [ i J + 2 . 0 * dyn ( i j ) ; 
f ree_vector (yt, 1, n) ; 
f ree^vector (dyt ,l,n); 
f ree_vector (dym, 1 , n) ; 
free_vector (dyn, 1, n) ; 

) 

/«*««****•****»«»**.•** rkdumb function ************************* 1 / 

void rkdumb (vs tart , r.var, xl,x2,nstep, derivs) 
int nvar,nstep; 
float vstart ( ] , xl , x2 ; 
void (*derivs ) ( ) ; 

{ 

int i, k; 
float x, h; 

float * v, *vout, *dv, *vectorr(); 
void rk4(), nrerror () , f r ee_vector ( ) ; 

v - vectorr (1 , nvar ) ; 
vout - vectorr ( 1 , nvar ) ; 
dv - vectorr (1, nvar ) ; 
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for (i-1 ; i<-nvar ; i++) 

( 

v|i] - vstart I i ) ; 
y U) U) - v | i ] ; 

) 

xx ( 1 ] - xl; 
x - x 1 ; 

h - (x2-xl ) /nstep; 
for (k-1 ; k<- nstep ; k++ ) 

{ 

( *derivs) (x,v,dv) ; 

rk4 (v,dv, nvar, x, h, vout , derivs) ; 

if (x+h — x) nrerror ( "Step sire too small in routine* RKnrnn") 

x ♦- h; 

xx[k+l] - x: 

for (i-1 ; i<-nvar ; i++) 

( 

v (i) - vout (i] ; 
y l i ] (k+ 1) - v(i); 

) 

) 

free_vector (dv, 1, nvar) ; 
free^vector (vout, 1 , nvar ) ; 
free^vector (v, 1, nvar) ; 

) 

/♦****«*************•* r kqc function *************************** 4 / 
void rkqc (y , dydx, n, x,htry,eps,yscal, hdid, hnext , derivs ) 

float ylJ, dydx ( ) , *x, htry, eps, yscal(], *hdid, *hnext; 
void ( *derivs ) ( ) ; 
int n; 

( 

int i; 

float xsav, hh, h, temp, errmax; 

float *dysav, *ysav, *ytemp, *vectorr () ; 

void rk4(), nrerror (), f ree_vector ( ) ; 

dysav - vectorr (1 , n) ; 
ysav - vectorr ( 1, n) ; 
ytemp - vectorr ( 1 , n) ; 
x s a v - ( * x ) ; 

for (i-1 ; i<-n; i + +) 

( 

ysav ( i ] - y [i] ; 
dysav ( i ) - dydx I i ) ; 

] 

h - h t r y ; 
for ( ; ; ) 

( 

hh - 0 . 5 * h ; 

rk4 (ysav, dysav, n,xsav,hh, ytemp, derivs) ; 

*x - xsav+hh; 

(♦derivs) (* x, ytemp, dydx) ; 

rk4 (ytemp, dydx, n, *x, hh, y, derivs) ; 

*x - xsav+h; 

if ( * x — xsav) nrerror ("Step size too small in routine > 

rk4 (ysav, dysav, n, xsav, h, ytemp, derivs ) ; 

errmax - 0.0; 

for (i-1 ; i<-n; i++) 

{ 

ytemp ( i] - y I i J -ytemp (i J ; 

temp - fabs (ytemp (i) /yscal ( i ) ) ; 

if (errmax < temp) errmax - temp; 
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) 

errmax /- eps; 
i f (errmax <- 1.0) 

I 

•hdid - h; 

*hnext - (errmax > ERRCON ? SAFETY *h * oxp ( PGPOW* 1 (o t i in. i >: ) ) 

4 n * ’-I ) 

break; 

) 

h - SAFETY*h*exp (PSHRNK*log (errmax) ) ; 

) 

for (i-1 ; i<-n; i++) y ( i ) +- ytemp ( i ] •FCOR; 
f ree_vector (ytemp, 1 , n) ; 
f ree_vector (dysav, i , n ) ; 
f ree_ vector (ysav, 1, n) ; 

) _ 

/••**«*******«***••** odeint function ***************** 4 *** / 
void odeint (ystart, nvar, xl,x2,eps,hl, hmin, nok, nbad, derivs, rkq<~) 

float ystartl], xl, x2, eps, hi, hmin; 
int nvar, *nok, *nbad; 
void ( *derivs ) ( ) ; 
void ( *rkqc ) ( ) ; 

l 

int nstp, i; 

float xsav, x, hnext, hdid, h; 
float *yscal, *y, *aydx, ‘vectorr () ; 
void nrerrorO, f ree_ vector ( ) ; 

yscal - vectorr (1 , nvar ) ; 
y - vectorr (1, nvar) : 
dydx - vectorr ( 1 , nvar ) ; 
x - x 1 ; 

h - (x2 > xl) ? f abs (hi ) : -fabs(hl); 

*nok - (*nbad) - kount - 0; 

for (i-1 ; i<-nvar ; i+*-) y ( i ] - ystart(i); 

if (kmax > 0) xsav - x-dxsav*2.0; 

for (nstp-1 ; nstp<-XAXSTP; ns tp++) 

{ 

('•derivs) (x,y,dydx); 

for (i-1 ; i<-r.var ; i++) yscal(i] - f abs (y ( i ) ) + f abs (dydx { i ) * h ) • 7 1 : J V ; 
if (kmax > C) 

( 

if ( f abs (x-xsav) > fabs (dxsav) ) 

l 

if (kount < kmax-1) 

{ 

xp[++kount) - x; 

for (i-1 ; i<-nvar ; i + + ) yp { i ] | kount ) -y i i ) ; 
xsav - x; 

) 

) 

) 

if { (x+h-x2) * (x+h-xl) > 0.0) h - x2-x; 

( *rkqc) (y, dydx, nvar, tx, h, eps, yscal, £ hdid, fihnext , derivs) ; 
if (hdid — h) ++{*nok); else ++(*nbad); 
if ( ( x-x2 ) * (x2-xl ) >- 0.0) 
l 

for ! i-1 ; i<-nvar; i+ + ) 

ystart (i) - y (i) ; 
if (kmax) 

( 

xp [ ++kount ) -x; 

for (i-1 ; i<-nvar ; i++) yp(i] (kount) « y(i); 
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) 

f ree_vector (dydx, 1, nvar) ; 
f ree_vector (y, 1, nvar ) ; 
f ree_vector (yscal f 1 , nvar ) ; 
return; 

} 

if (fabs(hnext) <- hmin) nrer ror ( "Step size too small in Onr.iMi ’*); 
h - hnext; 

) 

nrerror ("Too many steps in routine ODEINT" ) ; 

} 
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/ 



declarations 



♦include "rkk.h” 



♦define N 7 

♦define MAXARRAY 100000 
♦ define SQ(x) ( (x) * (x) ) 



int kount, kmax; 
float **yp, *xp; 
float dxsav, dxmin, dxlimit; 

float 11-0.0, 12-0.0, 13-0.0, M-0.0, m-0.0, mb-0.0, a=0.0, fc-0.0, 
float Ilp-0 .0, l2p-0 .0, 13p-0 . 0, Mp-0 . 0, mp-0 . 0, mbp-0 . 0, ap-0 . 0, kp-0 . 0. <*p 0 . «) ; 
float L-0.0, sigma-0. 0, HT-0.0, v-0.0, rho-0.0, rhop^O.O, L1-0.U, 1. 2-0.0: 

float A-0.0, C-0.0, J3p-0.0, w3p-0.0; 

/************************ EQUATIONS OF MOTION *************** ****.' 

void der ivs ( x, y , dydx) 
float x, y ( j , dydx [ ) ; 

{ 



float 


zeta-0 . 0, 


dydxzeta 


i-O.O; 


float 


Al-0.0, 


A2-0.0, 


A3-0.0, 


float 


A9-0.0, 


A10-0.0, 


All-0.0, 


float 


A17-0.0, 


A18-0.0, 


A19-0.0, 


float 


A25-0.0, 


A26-0.0, 


A27-0.0, 


float 


Bl-0.0, 


B2-0.0, 


B3-0.0, 


float 


B9-0.0, 


B10-0.0, 


Bll-0.0, 


float 


B17-0 .0, 


B18-0.0, 


B19-0.0, 


float 


B25-0.0, 


B26-0.0; 




float 


Cl-0.0, 


C2-C.0, 


C3-0.0, 


float 


C9-0.0, 


C10-0. 0, 


Cll-0.0; 


float 


Dl-0.0, 


D2-0.0, 


D3-0.0, 


float 


El-0.0, 


E2-C.0, 


E3-0.0, 


float 


E9-0.0, 


El 0-0 . 0; 




float 


Zl-0.0, 


Z2-0.0, 


Z3-0.0, 



A4-0.0, A5-0.0, A6-0 . 0, A7-0.0, 
A12-0.0, A13-0.0, A14-0.0, A1S-0.0, 
A20-0.0, A21-0.0, A22-0.0, A23-0.0, 
A28-0.0; 

B4-0.0, B5-0.0, B6-0.0, R7-0.0, 



B12-0.0, 

B20-0.0, 


B13-0.0, 

B21-0.0, 


B14-0.0, 

B22-0.0, 


B1S-0 .0, 
B23-O.Q, 


C4-0.0, 


C5-0.0, 


C6-0 .0, 


C7-0.0, 


D4-0.0, 


D5-0.0, 


D6-0.0, 


D7-0.0, 


E4-0.0, 


E5-0.0, 


E6-0.0, 


E7-0 . 0 , 


Z4-0.0, 


Fl-0.0, 


F2-0.0; 





Art-0. 0 

Al^-n.o 
A24-0 . n 



P3-0.0; 
niO-O . 0; 
1C ; 0.0; 



C‘» - 0 . O; 



n s - o. . o ; 

F* r . - 0 . O; 



dydx [ 4 ] - y ( 5 ] ; 



dydx [ 6] - y ( 7 ) ; 



zeta - rho*y [ 4 ] +rhop*y [ 6 J ; 



dydxzeta - rho* y ( 5 ] •» i hop * y 1 7 ] 



A1 

A3 

A5 

A7 

A9 

All 

A13 

A15 

A17 

A19 

A21 

A23 

A25 

A26 

A28 



- (A-C) *y[23*y(3); 


A2 


- 


J3p*sigma * y ( 2 J ; 


2 *MT* zet a ‘dydxzeta *y ( 1 3 ; 


A4 


- 


MT*SQ (zeta ) ; 


-MT*SQ ( zeta ) * y [ 2 ) *y ( 3 ) ; 


A6 


- 


-2*m* ( zeta«-L2 ) * y [ 4 ] ; 


2 *n* ( zeta + L2) *y ( 4 ) *y ( 2 ) *y [ 3 J ; 


A8 


- 


-2 *m* (zeta+L2) ^ y [ 5 ] * y [3 ] : 


-2*m*y ( 4 J *dydxreta*y ( 1) ; 


A10 


- 


2 *m*y (4 ) *y!5)*y(l]; 


m*SQ (y ( 4 3 ) ; 


A12 


- 


-m * SQ ( y [ 4 ) ) * y ( 2 ] * y | 3 } ; 


-m*y [ 4 ) *a; 


Al 4 


- 


-m*y (43*a*y(l]*y[2); 


- 2 * mp * (zeta-Ll) *y (6) ; 


Al 6 


- 


2 *mp* (zeta-Ll) *y [6 J *y [?) *y[3] 


-2 *mp* (zeta-Ll) *y [7 ) *y (1) ; 


Al 8 


- 


-2*mp*y ( 6 ) ‘dydxzola ' y ( 1 ] : 


2 *mp*y [ 6) *y ( 7 ) *y ( 1 ) ; 


A20 


- 


mp*SQ (y l 6) ) ; 


-mp‘SQ (y[6])*y[2)*y(3); 


A22 


- 


-mp * y ( 6 ) * ap * cos ( s i gm.n 4 >: ) ; 



-ir.p*y [ 6) *ap*cos (sigma*x) *y [ 1) *y [ 2 ] ; A24 - mp* ap * r, i :« ( s i qma 4 ) . 

mp*ap* sin (sigma *x) *y { 6J * (SQ (y [3) ♦ sigma) -SQ(y ( 2) ) ) ; 

A+A4 +A6+A1 1+A15+A20 ; A27 - A13+A22; 

A1*A2 + A3*A5 + A7+a8+A9 + A10 + A12+A14+A16 + A17+A18 + A^ A21-* A2 3» A? 5; 



B1 • - (C-A) *y ( l] *y ( 3] ; B2 - -J3p* sigma * y ( 1 ] ; 

B3 - 2*MT*zeta*dydxzeta*y (23 ; B4 - MT*SQ(zela); 
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B5 - MT*SQ(zeta) *y (1) *y(3) ; 

B7 - -2 *m* ( zeta+L2 ) *y [ 4 ] *y 1 1 ) * y 1 3 ] ; 
B9 - -2*m*y 1 4 ) *dydxzeta*y [2] ; 

Bll - m‘SQ(y (4) ) ; 

B13 m*a; 

B15 - -2*mp* (zeta-Ll) *y (6) ; 

B17 - -2*mp* (zeta-Ll) *y (7) *y (2) ; 

B19 - mp*SQ (y I 6) ) ; 

B2l ■» -mp*ap*cos (sigma*x) ; 



B6 - -2*m* ( zeta ' L2) *y |4 ; 

B8 • -2*m* ( zeta« L2) * y I 5 j • y I 2 | ; 

BIO - 2 *m*y I 4 ) *y I 5] *y I 7 | ; 

B12 - m*SQ (y I 4 ] ) *y | 1 ] ‘ y I 1 : ; 

B14 m*a*y ( 4 ) * (SO (y I 3 I ) - r O(y 1 1 II ) 

B16 - -2*mp* ( zeta-L) ) ‘ y ( 6 ‘ y I 1 I ' y I ’• I 
B18 - 2*mp*y (6) *y |2) * (y | 7 ; -<lyH:<:-«M .i) 
B20 - mp*SQ(y(6])*yt])*yP): 

B23 - -mp*ap* si n (siom * X )• y ((■>); 



B22 - -mp*ap*cos (sigma* x) * (SQ (y (3) +sigma) -SQ (y ( 1) ) ) *y ( 6) ; 

B24 - mp*ap*sin (sigma*x) *y ( 6) *y [ 1 j *y (2 ) ; B25 - A i B4 • HO ' : ) 1 < IU ' l> I" 

B26 - B1+B2+B3+B5+B7+B8+B9+B104B12+B14+B16+B17+B18+B20+B22+B24; 



Cl - -2*m*a*y (5) *y(l ) ; C2 - -m*a*yMJ; 

C3 - m*a*y (4 ) *y (2) *y (3) ; C4 - -2*mp*ap* sin (sigma *:•:)* y I 7 | • y |? 1 . 

C5- - -mp*ap*sin ( sigma*x) *y [ 6] ; C6 « -mp*ap*sin (sigma*x ) • y I f . ) * y I 1 I ‘ y M I : 

C7 - -2*mp*ap*cos (sigma*x) *y (7 ) *y ( 1 ) ; C8 ” -mp*ap*cos (sigma*x) • y I 6 1 ; 

C9 - mp*ap*cos (sigma*x) *y ( 6) *y (2 ) *y (3) ; 

CIO - C2+C8; Cll - C1+C3+C4 +C6+C7 4 C9 ; 



D1 - m*(l-rho); D2 - -mp*rho; 

D3 - -m* a ; D4 - m*a*y ( 1 ) *y (3) ; 

D5 - -m* (SQ(y(l) )+SQ(y (2) ) ) * (y (4) » (1-rho) -L2-rhop* y [6] ) ; P r -r*y|:.|; 

D7 - k*y[4); D8 - D4+D5+D6+D7; 

El - -m*rhop; E2 - mp* (1-rhop) : 

E3 - mp*ap*sin (sigma*x) ; E4 - mp*ap* sin ( sigma* x) *y ( 2 ) * (y I 3 ) 1 7 * ~ i gm i ) 

E5 - -mp*ap*cos (sigma*x) ; E6 - mp*ap*cos (sigma*x) * y ( 1 ) * ( y ! 3 ] ' 7 * r. i gir i ) 

E7 - -mp* (SQ(y(l) )+SQ(y [2] ) )* (y ( 6) * (1-rhop) +Ll-rho*y (4) ) ; F<? ^ *-p * y I 'l 

E9 - kp*y ( 6) ; E10 - E4+E6+E7+E0+E9; 

Z1 - (- (A27*B13) / (A26*B23) -E1/E3)/ ( (A27*B25) / (A26*B23) +E5/E3) ; 

Z2 - (-(A27*B21)/(A26*B23) +A24 /A26-E2/E3) / ((A27*B25)/(A26*B23) *F5/F.3) ; 

FI - (-(A27*B26)/(A26*B23) +A28/A26-E10/E3) / ( (A27*B25) / (A26*R23) * 2 5/E3) ; 

Z3 - ( (C*B13)/ (C10*B23)-(A27*B13)/ (A26*B23) ) / 

( (A27*B25) / (A26*B23)-(C*B25) / (00*073) *C5/ 010) : 
Z 4 - ( <C*B21)/ (C10*B23)- (A27*B21)/ (A26*B23 ) +A24/A2 6 ) / 

( (A27*B25)/(A26*B23)- (C‘B25) / (Cl 0 *B7 3 ) * CD/C) C ) ; 
F2 - ( (C*B26)/(C10*B23)-C11/C10-(A27*B26)/(A26*B23)+A28/A26)/ 

( (A27*B25) / (A26*B23)-(C*B25) / (CIO* L’ 73) <C5/C10); 



dydx [5] - ( (-F2-D8/D3)/ (Z4+D2/D3) + (F1+D8/D3) / (Z2+D2/D3) )/ 

( (-Z1-D1/D3)/(Z2+D2/D3H(Z3^1/P3)/ C4 11)7/1)3) ) : 



dydx (7) - ( (-F2-D8/D3)/ (Z3+D1/D3)+ (F1+D8/D3) / (Z1+D1/D3) )/ 

( (-Z2-D2/D3)/ (Z1 + D1/D3) + (Z4 < D2/D3) / ( 3 < 1)1 / 1 ' 7 ) ) ; 

dydx (3) - (dydx (5) * (C5*B13) / (C*B25) + 

dydx 1 7) * ( (C10*A24) / (C*A26 ) + {C5*B21 ) / (C* B2 5) ) + 

(C10*A28) / (C*A26) + (C5*B26 ) / (C*B25) -Cl 1/C) / 

(1- (C10*A27)/(C*A26)-(C5*B23)/(C*B25) ) ; 



dydx(l] - dydx ( 3) * (-A27/A26) +dydx (7) * (-A24/A26) + (-A28/A26) ; 

dydx (2) - dydx (3) * (-323/B25) +dydx (5) * (-B13/B25) +dydx ( 7 ] * ( -B21 /D75 ) -R2 6/ n2f, ; 

) 



MAIN PROGRAM / 

main ( ) 
i 

int i, j, nbad, nok; 

int 11-1000, 12-20C0, 13-3000, 14-4000, 15-5000, 16-6000, 17-7000; 
int 18-8000, 19-9000, 110-10000; 

float eps, hi, hmin, dxsavd, xl, x2, xscale, ‘ystart; 

float Hl-0.0, H2-0.0, H3-0.0, Hlp-0.0, H2p-0.0, H3p-0.0, zcm-0.0, n : gum'-iO ; 
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float Tl-0.0, T2-0.0, T3-0.0, T4-0.0, T5-0.0, T6-0.0, T/-0.0, 7°-'\0; 
float T9-0.0, T10-0.0, Tll-0.0, T12-0.0, T13-0.0; 

float zcmrigid, It, Itp, Ittotal, Is, Isp, Istotal, wn, lnmbdi. i-iTl-hp; 
float *ke, *kedel, *etotal, *etotaldel, 4 h, *hdel, *thota; 

float *Ec, *Ecdel, *Ecp, *Ecpdel, *Echype, *Echypedel, *Ecphyp~, * Frphypr** I** 1 
float Q, Qp, *eta, *etap, Ecmin, ecfactor, ecpfactor, hsq, Ecfiml, K.-pl iml 
float signec, signecp, signq, signqp, arg, argp, Ecf, F.cfp; 
float *Ecexp, *Ecpexp, *Ecexpdel, *Ecpexpdel, *etah, *otahp; 



ystart 


- 


vectorr (1,N) ; 


Ec 


- 


vector r (1 , y 


: AX APR AY ) 


ke 


- 


vectorr (1, MAXARRAY) 


Ecpdel 


- 


vector r ( 1 , MAX AP.P.AY ) 


kedel 


- 


vectorr (1, MAXARRAY) 


Echype 


- 


vectorr (1,1 


AX ARRAY ) 


etotal 


- 


vectorr (1, MAXARRAY) 


Echypedel 


- 


vectorr ( 1 , ? 


AX APR AY) 


etotaldel 


- 


vectorr (1, MAXARRAY) 


Ecphype 


- 


vectorr ( 1 , : 


ARRAY) 


h 




vectorr (1, MAXARRAY) 


Ecphypedel 


*■ 


vectorr ( 1 , 1 


,\,\:\RR A I ) 


theta 


- 


vectorr ( 1 , MAXARRAY) 


eta 


- 


vectorr ( i , : 


AX A PRAY ) 


hdel 


- 


vectorr (1, MAX ARRAY) 


etap 


- 


vector r ( 1 , ? 


A X ARRAY ) 


Ecexp 


- vectorr ( 1, MAXARRAY) ; Ecexpdel 


- 


vector r ( 1 , 1 


AX ARRAY) 


Ecpexp 


- vectorr (1, MAXARRAY) 


i ; Ecpexpdel 

etah 
etahp 


- 


vectorr ( 1 , : 
vectorr ( 1 , : 
vector r ( 1 , : 


AX ARP AY ) 
AX APR. AY ) 
A.X ARRAY ) 



scanf r%f %f %f%f%f%f%f%f%f%f%f%f%f%fmf%m%nf%muuuu?.f i:v. 

4x2,4ll,4l3,4M,4m,4a,*k,4c,4L,4llp,&l3p,4Mp,4mp,4ap,&kp. f . cp, 
tw3p, fcystart(l) ,4ystart[2],£ystart(3] , fiystart ( 4 ) , 4ysta rt (5J , 
4ystart[6],4ystart[7] , 4eps, 4 ecf actor , 4ecpf actor , 4 Ecf , 4E'*fp) : 

xl-0 .0001; hl-0. 00 001; hmin-1 . 0e-l 1 ; kmax-100000; 

mb - m; mbp - mp; 12 - II; I2p - lip; 

sigma - w3p-ystart ( 3) ; 

/**** critical damping - input damping - 1000.0, makes it crit: '■*»! *«»*/ 
if (c -- 1000. 0J c - 1 . 5* sqrt ( 4 *n*k ) ; 
if (cp — 1000.0) Cp - 1 . 5* sqr t ( 4 *mp* kp) ; 

/*** establish time interval for saving data ***/ 
dxmin - 1.0e-4; 
dxlimit - (x2-xl) /1750.0; 
dxsav - (x2<20.0) ? dxmin : dxlimit; 



MT - M+Mp+4 . 0 *mb + 4 . 0*mbp; 
v - (Mp+4 . 0*mbp) /MT; 

A - Il+Ilp+2 . 0*mb*SQ (a) +2 . 0*mbp*SQ (ap) * (Mp+4 .0*mbp) * (1 . 0-v) • r; C.) ; 
C - l3+I3p+4 . 0*mb*SQ (a) +4 . 0*mbp*SQ (ap) ; 

J3p - l3p+4 . 0*mbp*SQ (ap) ; 
rho - m/MT; 
rhop - mp/MT; 

LI - L* (M + 4 . 0 *mb) /MT; 

L2 - L* (Mp + 4 .0*mbo) /MT; 



zcmrigid 

It 



Itp 



Ittotal 

Is 

Isp 

Istotal 



( (Mp+ 4 . 0 *mbp ) * L) /MT; 

1 1 +M* SQ ( zcmrigid) +mb* (4 . 0*SQ ( 2 cmrigid) +2 . 0*SQ(a) ) ; 

I lp+Mp*SQ(L- zcmrigid) +mbp* (4 . 0*SQ(L- zcmrigid) *-2 . 0 * SC ) : 
It+Itp; 

13 + 4 . 0 *mb* SQ (a ) ; 
l3p+4 . 0*mbp*$Q (ap) ; 

Is+Isp; 



odeint (ystart,N,xl, x2, eps , hi , hmin, 4 nok , & nbad, der i vs , r kqc ) ; 

for ( j-1 ; j<-kount ; j++) 

l 

zcm - ( (Mp+4 . 0*mbp) *L+m*yp [ 4] ( j] +mp*yp[ 6) ( j ] ) /MT; 

HI - ( I1+M*$Q (zcm) +m* ( 2 . 0* SQ (a ) +4 . 0* SQ ( zcm) -t SQ (yp [ 4 ] ( j ; ■ 
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-2.0*zcm*yp[4) l j] ) ) *yp 1 1) I j) -m*a*yp I 4 I I j] *yp I 3 ] I j ) ; 

H2 - (I2+M*SQ (zcm) +m* (2.0*SQ(a)+4.0*SQ (zcm) +SQ (ypl 4 ] | j] ) 

-2 . 0* zcm*yp [ 4 ] I j] ) ) *yp|2] | j ] -m*a*yp I 5 ) I j) ; 

H3 - -m*a*yp[4) [ j] *yp(l) [ 3) + (13+4 . 0*m*SQ (a) ) *yp[3) I j] ; 

Hip - (Ilp+Mp*SQ (L-zcm) +mp* (2 . 0*SQ (ap) +4 . 0*SQ (L-zcm) + SQ (yp I 6] I j) H 
2 .0* (L-zcm) *yp|6) ( j] ) ) *yp|l) t j) -mp*ap*cos (sigma*xp| 3) ) * 
yp [ 6] [ j) *yp(3] I jj -mp*ap* sigma* cos (sigma*xp( j) ) * 

' yp 1 6) [ j)+mp*ap*sin(sigma*xpl j] ) *yp|7) I j) ; 

H2p - (I2p+Mp*SQ (L-zcm) +mp* (2 . 0*SQ (ap) +4 . 0*SQ (L-zcm) +SQ (yp I 6] | j) ) • 
2.0* (L-zcm) *yp| 6) I j) ) ) *yp[2] ( j)-mp*ap*sin(sigma*xp| j) ) * 
yp 16) ( j) *yp [3] ( j) -mp*ap* sigma* sin (sigma *xp ( j) ) * 
yp 1 6 ] 1 j ]-mp*ap*cos( sigma* xp[ j) )*yp(7) l j); 

H3p - -mp*ap*cos (sigma*xp[ j) ) *yp(6) | j) *ypll) ( J)-mp*ap* 

sin ( sigma* xp (j))*yp(6) (j)*yp(2) (j)+ (I3p+4 . 0*mp*SO (ap) ) * 

(VP 13 ) [ j)+sigma); 

h ( j) - sqrt (SQ (Hl+Hlp) +SQ (H2+H2p) +SQ(H3+H3p) ) ; 
argument - (H3+H3p) /h ( j) ; 
if (argument>0 . 99999999) argument « 1.0; 
theta(j) - 57 ,2957795131*acos (argument) ; 
hdel [ J) - h(3)-hll); 

T1 - 0.5*(Il*SQ(ypll)I3))+I2*SQ(ypt2] (3))+I3*SQ(yp(3) (3))); 

T2 - 0.5*(Ilp*SQ(yp(l) l31)+I2p*SQ(yp|2) I3))+I3p* 

SQ(yp|3] | j) * r. i 'jrm ) ) 

T3 - 0 . 5*m* (SQ(yp(5) ( 3 1 ) + (SQ (yp 1 1) ( 3 ] ) +SQ (yp [2 ) [ j] ) ) *SQ (yp M ) I i 1 ) 

-2 . 0*a*yp ( 1 ] I3)*yp|3] ( 3) ‘ypM 1 I j)) 

T4 - 0 . 5*mp* (SQ(yp[7] ( j))+SQ(yp[6) (j)) * (SQ(ypIl] [3))+SQ(yp|2] ( j!) ) 

+ 2 .0*ap*sin ( sigma *xp[ 3) ) * yp ( 2 ) l 3) *yp|3] [ 3 ] *ypl 6 ) [ j ) 
-2 . 0*ap*cos (sigma*xp[ 3 j ) * yp ( 1 J [ 3) *yp [3] ( 3) *ypl f >) I jl) 

T5 - -m*a*yp|2) (3)*yp[5) 13); 

T6 - -mp* (-ap*yp [7) | 3) *yp[2) | 3) *cos (sigma*xp [ 3) ) 

* ap*yp [ 7 ) I 3 J * yp 11) I3)*sin(sigma*xp(3)) 

-ap* sigma* yp [ 1 ) I j) *yp ( 61 (31 *cos (sigma*xp(3) ) 

-ap*sigma*yp (2) I j) * yp l 6) [ 3 ] * sin (sigma*xp| 3) ) ) ; 

T7 - (0.5*SQ(m*yp[5) (3)+mp*yp[7) (3) ) )/MT; 

T8 - (0.5* (SQ(yp[l) (3))*SQ(ypI2) (3)) ) 

*SQ (m*yp 1 4 ) [ 3 ) -<mp* yp | 6 ] [ j] )) /MT 

T9 - 0.5* (SQ(yp(l) 13) )+SQ(yp[2) (3) ) ) * (M+4 .0*m) *SQ(L) * 

SQ ( (Hp+4 . 0* mp) /MT) 

T10 - 0.5»(SQ(yp(l) I3))+SQ(yp|2) (3)))«(Mp+4.0*mp)«SQ(L)* 

SQ ( (M+4 . 0 *m) /MT) 

Til - (m+mp) * (-yp(5) 13) ) * (w*yp[5) (3)+mp*yp[7) [ 3) ) /MT; 

T12 - -m*yp|4) (3)*(SQ(yp[l) (3))+SQ(yp|2] 13)))* 

(m*yp [ 4 ) [ 3) + mp*yp (6) (3) + (Mp+4 . 0*mp) * I.) /MT 
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T 13 - -mp*yp(6] [ jlMSQ(yp(l] l j ] ) *SQ ( yp [ 2 ) [ j]) ) * 

(m*ypM J [ j ] “»mp*yp( f ] | 3 ] - (f. M 0 *m) 4 I.) /m 

Ke ( j j - T1+T2+T3+T4+T5+T6*T7*T8«T9+T10+T11*T1? <T1 i; 

kedellj] - ke(j]-ke(l); 

etotal ( j] - ke ( j] + 0 . 5 *k*SQ (yp( 4 ) [ j ] ) +0 . 5 *kp*SQ (ypl f>) l 1 1 ) ; 

etotaldel [ j ] - etotal [ j ] -etotal [ 1 ] ; 

Eclj] - 0 . 5 * (IttotalMSO(ypll] l j ] ) +SQ(ypl 2 ] [ j])) 

+Istotal*S 0 (ypl T] l i 1 ) ) 



Ecdellj] - Eclj]-Ecll]; 

Ecpl j] - 0 . 5 * (Ittotal* (S 0 (yplD l jl ) +SQ(yp[ 2 ) [ j] ) ) 

+ Is total*. SO (ypl 3 ] ( i) 'ninnn) ) 



Ecpdel [ j ] - Ecpl jJ-Ecpll] ; 

) 



wn - (Is*yp 13] 1 1 ] + lsp* (yp [3] U3 +sigma) ) /Ittotal; 

lambda - wn-yp( 3 ](l]; 

lambdap - wn- (yp [ 3] ( 1 ] +sigma) ; 

hsq - SQ(Ittotal) * (SQ(ypll) ll])+SQ(ypl 2 ] 11 ] ) ) 

+SQ ( Is *yp 13 ] HJ+Isp* (yp [ 3 ] { 1 ] minim) ) 



if (theta (kount ] -theta l 1] < 0.0) /**** stable condition *•**/ 

{ 

if (lambda < 0.0) { signec - 1.0; signq « 1.0; ) 

else l signec - -1.0; signq - -1.0; ) 



if (lambdap < 0.0) 
else 
1 

else 

{ 

if (lambda < 0.0) 
else 



( signecp - 1.0; signqp -1.0; ) 

{ signecp - -1.0; signqp « -1.0: ) 

/**** unstable condition ♦**♦/ 

( signec - -1.0; signq *1.0; ) 

( signec - 1.0; signq - -1.0; ] 



if (lambdap < 0.0) { signecp - -1.0; signqp *1.0: ) 
else ( signecp - 1.0; signqp * -1.0; ) 
1 



Ecfinal - (Ecf — 0 . 0 ) ? Ec (kount) : Ecf; 
Ecpfinal - (Ecfp — 0 . 0 ) ? Ecp(kount) : Ecfp; 



for ( j-1; j<-kount; j + + ) 

( 

Echypel j) - Ec(l] *Ecf inal / (Ec [ 1 ] + (Ecf inal-Ec ( 1 ) ) * 

exp (ecf actor * Ecf iml *:<p ( j ) ) ) 



Echypedel [ j ) - Echype [ j ] -Echype [ 1 ] ; 

Ecphype ( j ] - Ecp(l] *Ecp final / (Ecp [ 1 ) + (Ecpfinal -Ecp 1 1 )) * 

exp (ecpfactor *Ecpf mil 4 yp ( j ) ) ) 



Ecphypedel ( j ) - Ecphype ( j ] -Ecphype [ 1 ] ; 

Q - 2 . 0*Echype l j ) * Istotal* (1 . 0- ( Istotol/ Ittotal ) ) 

-hsq* (Istotal/Ittotal) * (1.0 - (Istotal/Ittotal) ) 
+(Istotal/Ittotal)*SQ(Isp) * SO (sigma) ; 

Qp - 2 . 0* Ecphype [ j) *Istotal* (1 .0- (Istotal/Ittotal ) ) 

-hsq* (Istotal/Ittotal) * (1 .0- (Istotal/ It total) ) 
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+ (Istotal/Ittotal) *SQ(Is) *SQ(niqmn) ; 


arg 


- ( ( Ittotal/ ( I ttotal -I slota ) ) ) 

* ( ( Isp*s i gma ) ♦ (s ignq 4 sqr 1. (0) ) ) ) / ( 1 1 

if (arg > 0.99999999) arg-1.0; 


etah l j) 


- 57 . 2957795131 *acos (arg) ; 


argp 


- (( Ittotal/ (Ittotal-Istotal ) ) 

* ( (-Is*sigma) « (signqp 4 sqrt (Op) ) ) ) / (- ;, r * Or p 1 
if (argp > 0.99999999) argp~1.0; 


etahp ( j ] 


- 57 . 2957795131 *acos (argp) ; 


Ecexp ( j J 


- (Ecll)-Ecf inal) * exp (ecf actor • xp ( j ) ) <-Rc f i n-i * : 


Ecexpdel l j) 


- Ecexp ( j ) -Ecexp l 1 ) ; 


Ecpexpl j) 


- (Ecpl 1) -Ecpf inal) *exp(ecp£actor* >:p | j) ) «F>pf inn 1 ; 


Ecpexpdel ( jj 


- Ecpexpl j) -Ecpexpl 1] ; 


Q 


- 2 . 0*Ecexpl j) *Istotal* ( 1.0- (Is total /I ttoi.nl ) ) 
-hsq* (Istotal/Ittotal) 4 (1.0- ( Istota 1/ I ttotn 1 ) ) 
+ (Istotal/Ittotal) *SQ(Isp) *SQ(sigma) ; 


Q? 


- 2 . 0* Ecpexpl j ) * Is total * (1.0- (Istotal/2 ttotn ) ) ) 
-hsq* (Istotal/Ittotal) 4 (1.0- ( Istotol/ II t^t n 1 ) ) 
+ (Istotal/Ittotal) *SQ(Is) *SQ( sigma) ; 


arg 


- ( (Ittotal/ ( Ittotal-Istotal) ) 

* ( (Isp*sigma) + ( s ignq * sqr t (0) ) ) ) / (r.qr i ) 

if (arg > 0.99999999) arg-1.0; 


cta( j) 


- 57 . 2957795131*acos (arg) ; 


argp 


- ( (Ittotal / (Ittotal-Istotal) ) 

* ( (-Is* sigma) + (signqp*sqrt (Qp) ) ) ) / (r.qr i (hsq) ) 
if (argp > 0.99999999) argp-1.0; 


ctapl j) 
) 


- 57 . 2957795131*acos (argp) ; 


/ * * * • 




init(); color_scale ' 
grey_scale( "greyscale 
bgcol (7) ; erase () ; 
nove (220, -110) ; print 


“cyanblue") ; 

1") ; windowO ( ) ; 
rolor (0) ; 

i ( "KINGORI DUAL SPINNER MASS SPRING SYSTEM") ; 



window (22, -85, 979, 872 ; bgcol(3); 



erase (); scale (0,10- 
size (1, 1) ; 


:3, 10000,0) ; rect (0, 0, 10000, 10000) ; 

rove (20, 110+200) ; date () ; 

r:ve (-160, 110+70) ; printf ("0") ; 

rove (10, 110-90 ) ; printf ( H w3 - w3 initial**); 


vector (0, 15, 110, 19) ; 


rove (-160, 15+70) ; printf ( "0") ; 
rove (10,19 + 200); printf ("w2") ; 


vector (0, 13, 110, 18) ; 


rove ( -160, 18+70) ; printf ("0") ; 



rove (10, 18+200) ; printf ("wl" 
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vector (0, 1 7 , 110, 17) 
vector (0,16,110,16) 
vector (0, 15,110,15) 
vector (0, 11,110,14) 
vector (0, 13,110,13) 



move (-160, 17470) ; 
move ( 10, 17 + 200) : 

move (-160, 16+70) ; 
move ( 10 , 16 ♦ 200) ; 

move (-160, 15 + 70) ; 
move ( 10, 15 + 200) ; 

move (-160, 14+70) ; 
move (10, 14 + 200) ; 

move (-160, 13 + 70) ; 
move (10, 13+200) ; 



print f ( "0") ; 
print £ ("z prime d 

printf ("0") ; 
printf ( "z prime") 

printf ("0") ; 
printf ("z dot") ; 

printf ("0") ; 
printf ("z") ; 

printf ("0") ; 
printf ( "ke-ke i 



) : 



( r ) , o - o \ ( i > ) ” ) : 



vector (0, 12,110,12); move (-160, 12+70); 

move ( 10, 12 + 200) ; printf ("Ec-Ec 



printf ("0") ; 
i (r), Echypo-Echype 



(h) ") : 



vector (0, 11,110,11) ; move (-160,11+70) ; printf ("0") ; 

move ( 10, 11 + 200) ; printf ("Ecp-Ecp i (r), Ecphype-Ecphypo i (•’)"); 



vector (0, 0, 110, 0) ; move (-160, 0+70) ; printf ("0"); 

move ( 10, 0+200) ; printf ( "theta, eta (b) , etap (r), h - h imi i.+l"): 



vector (0, 19+750, 11, 19+750) 
vector (0, 19+500,11, 19+500) 
vector (0, 19+250,11, 19+250) 
vector (0, 19-250, 11,19-250) 
vector (0, 19-500, 11,19-500) 
vector (0,18+250,11, 18+250) 
vector (0, 18-250, 11, 18-250) 
vector (0,18-500, 11, 18-500) 



vector (19,19+750, 110, 19+750) 
vector (19,19+500,110, 19+500) 
vector (19,19+250,110, 19+250) 
vector (19,19-250, 110, 19-250) 
vector (19, 19-500, 110,19-500) 
vector (19,18+250, 110, 18 + 250) 
vector (19,18-250, 110, 18-250) 
vector (19, 18-500, 110, 18-500) 



xscale - 10000 . 0/xp[kount ) ; 

color (1); /****• blue - total energy, rotor hype *****/ 

for ( j-1; j<kount; j++) 

( 

vector ( (int) (xscale*xp [ j] ) , (int) (13+ (0.01*etotaldel l j) ) ) , 

(int) (xscale* xp[ j + 1] ) , (int) (13+ (0 . 01*etotaldel ( j<* 1] ) ) ) ; 
vector ( (int) (xscale *xp( j) ) , (int) (13+ (0. l*etotaldel ( j) ) ) , 

(int) (xscale* xp ( j+1] ) , (int) (13+ (0.1 *etotaldel ( j+l ] ) ) ) : 



vector ((int) (xscale*xp(j]), (int) (12+(50 .0*Echypedel ( j ] ) ) , 

(int) (xscale* xplj + 1]), (int) (12+(50. 0*Echypcdcl ( j + 1 ) ) ) ) ; 

vector ( (int ) (xscale *xp ( j] ) , (int) (0+ (100 . 0*etah ( j ] ) ) , 

(int) (xscale*xp ( j+1 ] ) , (int) (0+ (100. 0*etah (j+1]))); 
vector ( (int) (xscale *xp ( j] ) , (int) (11+ ( 1000 . 0 *etah ( j ] ) ) , 

(int) (xscale*xp(j+l]J, (int) (ll+(1000.0*etah(j+l]))); 



1 



color (5); /***.* purple - platform exponential*****/ 

for ( j-1; jOcount ; j++) 

( 

vector ( (int) (xscale*xp(j)), (int) <11+ (50 . 0*Ecpexpdcl ( j] ) ) , 

(int) ( x scale *xp( j + 1 ) ) , (int) (11+ ( 50 . 0 *Ecpexpdel | j • ) j ) ) ) ; 

vector ( (int) (xscale *xp( j] ) , (int) (0+ (100 .0*etap( j] ) ) , 

( int ) (xscale*xp[j + l)), (int) (0+ ( 100 . 0*etop (j<l)))); 
vector((int) (xscale*xp(j]), (int) (ll+(1000.0*etap(j))), 

(int) (xscale*xp(j + l)), (int) (11+ (1000. 0*etap ( j « 1 ) ) ) ) : 
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) 



color(6); /**••• yellow - rotor exponential* 444 */ 

for ( j-1: J<kount; j + -) 

( 

vector ( (int) • xscale 4 xp ( j ) ) , (int) (12 + ( 50 . 0 4 Ecexpdc 1 I j) ) ) , 

(int) (xscale 4 xp ( j + 1 ) ) , (int) (12+ (50 . 0 * Ecexpdc 1 | i * t I i ' > 



vector ( (int ) [xscale 4 xp( j]) # (int) (0+(100.0*eta[j))), 

• (int ' (xscale 4 xp ( j + 1 ) ) , ( int ) (0 + ( 100 . 0 4 eta ( j-* 1 I ) ) ) ; 
vector ( (int; ixscale 4 xp( j] ) , (int) (ll+(1000.0*cta(j))), 

(int ' (xscale 4 xp ( j + 1 ) ) , (int) (11** (1000 . 0 4 eta | j < 1) ) ) ) : 

) 



color(0);* / 44 ** 4 black - platform hyperbolic * 44 * 4 / 

for ( j-1; J<kount; j+**I 

{ 

vector ( (int) ( xscale 4 xp ( j ]) , (int) (ll+(50. 0 4 Ecphypedo 1 ( j ) ) ) , 

(inti (x scale 4 xp (j + 1) ) , (int) (11+ (50 . 0* Ecphypcde l (j'l)))): 

vector ( (int ) ( xscale 4 xp( j] ) , (int) (0+(100.0*etahp(j))), 

(int ) (xscale 4 xp (j + 1)), (int) (0+(100.0*etahp(j+l)))); 
vector ( (int ) ( xscale 4 xp [j)), (int) (11 + (1000.0 4 etahp ( j ) ) ) , 

(int) (xscale 4 xp (j+1)), (int) (11+ (1000. 0 4 otahp( j '!)))): 



color (4); /***-• re d - rotor 44444 / 

for ( j-1; j<kount; j++) 

( 

vector ( (int) (xscale 4 xp ( j ) ) , (int) (110 +(14 32. 39448783* (yp(3) ( j) -yp!3J ( 1 ) ) ) ) 
(int) (xscale 4 xp( j+1) ) , (int) (110+ (1432 . 394 4 8783* (yp ( 3 ) ( j + 1) -ypl 3] ( 1 ) ) ) ) ) 

vector ( (int) (xscale 4 xp( j)), (int) (19+ (14 32. 39448783* yp (2) (j))), 

(int) (xscale 4 xp( j + 1)) , (int) (19+(1432.39140783 4 yp[2) ( j < 1 J ) )) ; 

vector ( (int) (xscale 4 xp( j)), (int) (18+(1432. 39448783* yp ( 1 ) (j)) ) , 

(int) (xscale 4 xp[ j + 1)) , (int) (18+ (1432 . 394 48783*yp ( 1 ) ! i H 1) ) ) : 

vector ( (int) (xscale 4 xp ( j) ) , (int) (17+ (1000 4 yp( 7) ( j ) ) ) , 

(int) ixscale 4 xp(j+l)), (int) (17+(1000 4 yp(7](j+l)))); 

vector ( (int) (xscale 4 xp ( j) ) , (int) (16+ ( 1000 4 yp ( 6) ( j ] ) ) , 

(int) [xscale 4 xp (j+1)), (int) (16+ (1000* yp (6) (j+1)))); 

v' or ( (int) (xscale 4 xp( j]), (int) (15+(1000 4 yp(5) (j))), 

(int) (xscale 4 xp (j + 1)), (int) (15+ ( 1000* yp ( 5 j ( j ♦ 1 ] ) ) ) ; 

vector ( (int) * xscale 4 xp( j)), (int) (14+(1000 4 yp(4) (j))), 

(int; [xscale 4 xp ( j + 1 ) ) , (int) (14+(1000* yp (4) (j'l)))): 

vector ( (int) ( xscale 4 xp ( j) ) , (int) (13+ (0 . 01 4 kedel ( j) ) ) , 

(int) xscale 4 xp ( j + 1 ) ) , (int) (13+ (0 .01 4 kcdel ( j + 1 ) ) ) ); 

vector ( (int) (xscale 4 xp ( j ) ) , (int) (13+ (0 . l 4 kedel ( j ) ) ) , 

(int) [xscale 4 xp ( j + 1 ) ) , (int) (13+ (0 . 1 4 kedel ( j*» 1 ) ) ) ) ; 

vector ( (int) ( xscale 4 xp( j) ) , (int) (12+ (50 . 0*Ecdel ( j ) ) ) , 

(int) .xscale 4 xp (j + 1)), (int) (12+ (50 . 0*Ecdel (j'l)))): 

vector ((int) (>:scale 4 xp ( j ) ) , (int) (11+ (50 .0*Ecpdel ( j ) ) ), 

(int) [xscale 4 xp( j+1)), (int) (11+ (50 . 0 4 Ecpdel ( j « 1 ) ) ) ); 

vector ( (int) (>:scale 4 xp( j) ) , (int ) (0+ (1000. 0*hdcl I j ) ) ) , 

(int) 'xscale 4 xp( j + 1)), (int) (0+0000 . 0*hdo) (j'l)))): 
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vector(tint) (xscale 4 xp|jj), ( int ) (0« ( 100 . 0 4 Lhct a | j J ) I , 

( int ) (xscale 4 xp(j*«l )) , (int) (0* ( 100 . 0* theta ( j • t 1 1 I ) •* 
vector ( (int) (xscale *xp ( j) ) , (int) ( 1 1 *• ( 1000 . 0 * thct «i I i ) ) ) , 

(int ) (xscale*xp(j+l)), (int) (11 * (1000 . 0M beta I j * 1 ) I I I ; 

vector((int) (xscale 4 xp ( j )) , (int) (0-»(0.0001 4 kc(j))), 

(int) (xscale 4 xp ( j+1 ] ) , (int) (0+(0.0001 4 kelj»l))) ); 

vector((int)(xscale 4 xp(j)), (int) (0+(0.0001 4 etotal|j))), 

(int) (xscale*xp( j + 1 ] ) , (int) (0* (0.0001 4 etotal | ) « I I ) ) ) : 

vector ( (int) (xscale*xp[ j] ) , (int) (0+ (0 . 1 *h I j ] ) ) , 

(int) (xscale* xplj+1)) , (int) (0*(0.1*h|j+l))) ); 



windowO ( ) ; 
color (0 ) ; 

/♦** print initial, final data 4 
move (350, 675) ; printf(" INITIAL CONDITIONS 



kount _ * cl" , I onn! ) 



move (22, 910) ; print f ( "time -%3.0f seconds”, xp 1 kount J ) ; 
move (500, 910) ; printf ("eps - \e %d/%d" , eps , nbad, nok) ; 

move (22,940); 

printf (**W1 - % 4 . 2 f W2 -\4.2f W3 4 . 2f " , yp [ 1 ) [ 1 ] , yp [ 2 ) ( 1 ) , yp I 3 ) I 1 ) ) : 
move (500, 940) ; pr int f ( "w3p- %9.7f L- %3 . 1 f M , w3p, L) ; 



move (22, 980) ; prir.tf("Il -%4.0f 12 -%4.0f 13 -%4 . Of M , II , 12, M) 

move (22,1010); printf ( N M - %4.0f m -%7 . 5f ", M, m) ; 

k -%4.1f c -%4 . 1 f " , a, k , c) ; 
zdot-%4 . 2f “, yp ( 4 ) ( 1 ], yp ( 5 )(!)); 



move (22, 1040) ; printf("a -\4.1f 
move (22 , 1070) ; printf ("z -\4.2f 



move (500, 960) ; print f ( "I lp-%4 . Of l2p-%4 .Of l3p-\4 .Of", Up, I2p, I 3 j>) : 
move (500, 1010) ; printf (”Mp -%4.0f mp-%7 . 5f ",Mp,mp) ; 

move ( 500, 104 0) ; print f ( "ap-%4 . If kp-%4.1f cp-*4 . If ", ap, kp, cp) ; 

move (500, 1070) ; printf (”zp -%4.2f zdotp -%4 .2f ", yp [6] ( 1) , yp[ 7) | 1 ) ) ; 

move (22,1110) ; printf( M h i-%9.4f ke i-%7. If", h(l] , ke [ 1] ) ; 
move (22 , 1135) ; printf(”h f-%9.4f ke f-%7 . If ", h (kount) , kc [ kount )) ; 
move (22, 1160) ; printf ("hdelper-%5 . 2f ", 100.0 4 (h ( kount ] -h | 1 ) ) /h| 1) ) ; 
move (250, 1160) ; print f ("kedelper-%5 . 2f", 100.0 4 (ke l kount ] -ke 1 1 ) ) /Y.* (1)1; 
move (500, 1160) ; 

print f (”e celper-%5 . 2f " , 100.0*(etotal(kount)-etotal(l))/^Lotn) |) )); 



move (730, 1160); 


printf ("£cp 


delper-%5 . 2f ", 100.0* (Ecp I kount ) -Erp 1 1 ] ) /Ecpl 1 I ) 


move ( 500 , 1110) ; 


printf ("Ec 


i-\9 


. 4 f theta i-%7 . 5f", Ec ( 1 ), them ( 11 ) ; 


move (500, 1135) ; 

*/ 


printf ("Ec 


f-%9 


. 4 f theta f-%7 . 5f", Ec ( kount ), thctn | konni ) ) 


move (22,900); 


printf (" 


platform rotor") : 


move (22 , 920 ) ; 


printf ("Up 


- 


% 6 . 2 f ”, Up) ; 


move (2 40, 920); 


printf ("11 


- 


%6.2£",I1); 


move (22, 940) ; 


printf ("I3p 


- 


\6.2t" , I3p) ; 


move (240, 940) ; 


printf ("13 


- 


% 6 . 2 f " , 13); 


move (22 , 9 60 ) ; 


printf ("Mp 


- 


l6.2£",Mp); 


move (2 4 0,960); 


printf (”M 


- 


%6 . 2f " , M) ; 


move (22, 980); 


printf ("mp 


- 


\6.2£”,mp) ; 


move (240, 980) ; 


printf ("m 


- 


%6 . 2f ", m) ; 


move (22, 1000); 


printf ("ap 


- 


1 6 . 2 £ " , ap ) ; 


move (240, 1000) ; 


printf ("a 


- 


% 6 . 2 f " , a ) ; 


move (22, 1020); 


printf ("kp 


- 


16 . 2f ", kp) ; 


move (240, 1020 ) ; 


printf ("k 


- 


%6.2£”,k) ; 
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move (22, 1040) ; 


print f ( "cp 


- 


\6. 2C", cp) ; 


move (240, 1040) ; 


printf ("c 


- 


% 6 . 2 f " , c ) : 


move ( 22 , 1060) ; 


print f ( "L 


- 


%6 . 4 f " , I,) ; 


move (240, 1060) ; 


printf ("Is/It 


— 


%6.3f", Istotal/Ittotnl) ; 


move (22,1100); 


printf ("zp 


- 


%6 . 2f ", yp I 6 J (3)); 


move (240, 1100) ; 


print f ( "z 


m 


% 6 . 2f ” , yp [4 J (I)); 


move (22, 1120) ; 


printf ( "dzp 


" 


%6 . 2f ", yp ( 7 ) |1J ) ; 


move (240, 1120) ; 


printf ( "dz 


- 


%6 . 2f ", yp [ 5) ID); 


move (22, 1140) ;• 


print f ( "w3p 


- 


% 6 . 2 £ " , yp I 3 J 1 1 J 4sigmn) : 


move (240, 1140) ; 


printf ( "w3 


- 


% 6 . 2f ", yp [ 3 J ID); 


move (22,1160); 


printf ( "wl 


- 


% 6 . 2f ", yp ( 1 ) ID) ; 


move (240, 1160) ; 


printf ( "w2 




%6.2f",yp|2J ID): 


move (480, 900) ; 


printf ("%3-0f 


seconds ", xp ( kount J ) ; 


move (730, 900 ) ; 


printf ("h i 




- % 9 . 4 f " , h [ 1 ) ) ; 


move (480, 920) ; 


printf("%d / 


%d 


/ %d" , nbad, nok, kount. ) ; 


move (730, 920) ; 


printf("h f 




- %9 . 4f'\ h [kount] ) ; 


move (480,940); 


printf ( "eps 




- \7e", eps) ; 


move (730, 940) ; 


printf ( "hdelp 




- %5.3f", 100.0* (h [ kount) -It 1 1 1 ) / i: i 1 ! 


move (4 80, 960) ; 


printf("ke i 




- %9 . 4£", ke 1 1 ) ) ; 


move (730, 960) ; 


printf("e i 




- %9 . 4 f " , etotal 1 1 ) ) ; 


move (480, 980) ; 


printf ("ke f 




- \9 . 4f " , ke [ kount ) ) ; 


move (730, 980) ; 


printf ( "e f 




- \9 . 4f", etotal I kount )) .- 


move (4 80, 10 00 ) ; 


printf ( "kedelp 


- %5. 3f", 100.0* (ke | kount ) -kc 1 ) ] ) / »•« 


move (739, 1000) ; 

'printf ( 


"edelp - %5. 


3f 1 


", 100.0* (etotal I kount I -etotal 11))/^' 


move (480,1030); 


printf ( "Ecpi 




- %9 . 3f ", Ecp( 1 ) ) ; 


move (730, 1030) ; 


printf ("Ec i 




- %9 . 3f ", Ec [ 1) ) ; 


move (480, 1050) ; 


printf ( "Ecpf 




- %9 . 3£" , Ecp ( kount )) ; 


move (730, 1050) ; 


printf ("Ec f 




- %9.3£", Ec (kount)) ; 


move ( 4 80, 1070 ) ; 


printf ("lambdap 


- % 5. 3£", lambdap) ; 



move (730, 1070) , 
move (4 80, 10 90) , 
move (730, 1090) , 
move (480, 1110) ( 
move (730, 1110) , 



print f ( "lambda 
print f ("signqp 
printf ("signq 



%5 . 3f ", lambda) ; 

%5 . 2f ", signqp) ; 

\5.2f ", signq) ; 
printf ( "d Ecp/ lambdap- \6. 3f ", (Ecp (kount ) -Ecp I 3 ))/ !.“*•: 
printf ("dEc/lambda - %6.3f", (Ec [ kount ) -Ec ( 1 ))/ 1 nirW 



move ( < 8 0 , 1130) ; print f ("Ecfp 
move(730, 1130) ; printf ("Ecf 



- %8 . 6f ", Ecfp) , 

- %8.6f",Ecf) ; 



move (480, 1155) ; 
move (620, 1155) ; 
move (760,1155); 
move (480,1175) ; 
move (620, 1175) ; 
move (760, 1175) ; 



print f ( "theta i -%6 . 2f", theta ( 1 )) ; 
printf("etap i -%6 . 2f ", etap [ 1 ) ) ; 
printf("eta i -%6 . 2f ", eta ( 1) ) ; 
print f ( "theta f -%6 . 2 f" , theta [ kount ]) ; 
printf ("etap f -% 6 . 2 f ", etap ( kount ) ) ; 
printf ("eta f -%6 . 2 f ", eta [ kount ) ) ; 



move ( 600, -33 ) ; 
move ( 600, -13) ; 
move ( 600, 27 ) ; 
move (600, 47) ; 



/ *** upper right hand corner 
printf ("Ec initial/final - \g 
printf ("Ecp initial/final - %g 
printf ("Ec factor / Ecf - ^g 
printf ("Ecp factor / Ecfp - \q 



***/ 

/ %g" , Ec | 1 ) , Ec I kount ] ) 
/ %g" , Ecp ( 1 J , Ecp [ V.eurU 
/ %g", ecf actor, Ecf} ; 

/ tg", ecpfactor , Ecf p) ; 



move (980, -33) ; printf ("%d" , nbad) ; 
move (980,-13) ; print f ( "%d" , nok ) ; 
move (980, 7) ; print f (" %d" , kount ) ; 
move (980, 27) ; printf ("nbad" ) ; 
move (980, 47) ; print f ( "nok " ) ; 
move (980, 67) ; print f ( "kount " ) ; 

move (980, 97) ; printf ("\6 . 4f ", (100 .0* (ke ( kount ] -ke l D ) ) /ko[ 1 ) ) ; 



Mi). 

1 1 | i 



rl'dnp) 
!-i ) ; 
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move (980, 117) ; pr int f ( " \ 6 . 4 f " , ( 1 00 . 0 * (Ctotnl I kount )-o total ! 1) >)/•’''» ^ M ) ) 
move (980, 137) ; pr int f C'\6 . 4 l " , ( 100 . 0 4 (h| kount. 1 -h | 1 )) ) /h | 1 ) ) ; 

move (980,157); printf ("kedel") ; 
move (980, 177) ; printf ("edelp") ; 
move (980,197); printf ("hdelp") ; 

move (980, 227) ; print f ("%g", theta (1) ) ; 
move (980, 247) ; printf ("%g", theta (kount ) ) ; 

move (980, 267) ; . printf ( M th 1") ; 
move (980, 287) ; printf ("th f ") ; 

move (98 0, 317) ; printf ("%g", eta ( 1] ) ; 
move (980,337); printf ("%g", eta | kount ) ) ; 

move (980, 357) ; printf ("eta 1"); 
move ( 980 , 377 ) ; printf ("eta £"); 

move (980, 407 ) ; printf ("%g", etap ( 1 ) ) ; 
move (980, 427) ; printf ("%g", etap (kount) ) ; 

move (980, 447) ; printf ("etapl") ; 
move (980,467); printf ("etapf ") ; 

move (980, 4 97) ; printf ("%g", Ec ( kount ) -Ec [ 1 )) ; 
move (98 0,517); printf ("tg", lambda) ; 
move (980, 537) ; printf ("%g", signq) ; 

move (98 0, 557) ; printf ("Ecdel") ; 
move (980,577); printf ("lambda") ; 
move (980, 597) ; printf ("signq") ; 

move (980, 627); printf ("%g", Ecp| kount) -Ecp ( 1 ) ) ; 
move (980,647); printf ("%g", lambdap) ; 
move (980, 667); pr int f ("%g", signqp) ; 

move (980, 687) ; printf ("Ecpdel") ; 
move (980,707); printf ("lambdap") ; 
move (980, 727); printf ("signqp" ) ; 

move ( 980, 757 ) ; printf ("%g", Istotal/Ittotal) ; 
move (980, 777) ; printf ("Is/It") ; 

f ree_matrix (yp, 1, N, 1, MAX ARRAY) ; 
free_vector (xp, 1 , MAXARRAY) ; 
f ree_vector (ystart,l,N) ; 
free_vector (ke, 1, MAXARRAY) ; 
free_vector (kedel, 1, MAXARRAY) ; 
free_vector (etotal, 1, MAXARRAY) ; 
free_vector (etotaldei, 1, MAXARRAY) ; 
free_vector (h, 1 , MAXARRAY) ; 
f ree_vect or (theta, 1, MAXARRAY) ; 
free^vector (hdel, 1, MAXARRAY ) ; 
f ree_vector (Ecexp, 1, MAXARRAY) ; 
free_vector (Ecpexp, 1, MAXARRAY) ; 
free_vector (etah, 1, MAXARRAY) ; 

) 



free_vector (Ec, 1, MAXAPRAY) ; 
free__vector (EcdcJ , 1 , MAX Ai r ') : 
free_vector (Ecp, 1, MAXARRAY ) ; 
f ree_vect or (Ecpdel , 1, MAXARRAY ) . 
free_vector (Echype, 1, MAXARRAY) : 
free_vector (Echypcdcl, 1 , MAX API* AY ) ; 
free_vector (Ecphype, 1 , MAX ARRAY ) ; 
f ree_vector (Ecphypcdel, 1 , MAXARRAY) ; 
free_vector (eta, 1, MAXARRAY) ; 
free_vector (etap, 1 , MAXARP.AY ) ; 
f ree_vector (EcexpdeJ , 1, MAXARRAY) ; 
free_veCtor (Ecpexpdol, 1 , MAXAKR.AY) ; 
free_vector (ctahp, 1, MAXARRAY) : 



132 



INITIAL DISTRIBUTION LIST 



1 . Defense Technical Information Center 2 

Cameron Station 

Alexandria, Virginia 22304-6145 

2. Library, Code 52 2 

Naval Postgraduate School 

Monterey, California 93943-5002 

3. LCDR V. M. Ortiz 2 

917 Copper Stone Circle 

Chesapeake, Virginia 23320 

4. Prof. I. M. Ross 8 

Code AA/Ro 

Department of Aeronautics and Astronautics 
Naval Postgraduate School 
Monterey, California 93943-5002 

5. Prof. H. A. Dahl 1 

Code PH/Dh 

Department of Physics 
Naval Postgraduate School 
Monterey, California 93943-5002 



133 



OUDLEY KNOX LIBRARY 
NAVAL POSTGRADUATE SCHOOL 
MONTEREY CA 8394S-5101 



